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ABSTRACT 


A  block  of  metal  is  subjected  to  a  concentrated  heat  source  resulting  in  a  pool  of 
molten  metal  surrounded  by  a  portion  of  the  unmelted  metal.  The  governing  system  of 
equations  is  known  as  the  weldpool  problem.  The  weldpool  problem  is  discretized  using  finite 
differences  to  discretize  time  derivatives  and  the  Finite  Volume  Element  method  (FVE)  to 
discretize  spatial  derivatives.  Multigrid  methods,  known  to  be  eflFective  on  uniform  grids, 
make  use  of  overlapping  uniform  grids  of  different  scales  to  better  approximate  solutions 
to  the  weldpool  problem.  However,  the  solid-liquid  interface  requires  extremely  fine  grid 
resolution  apd  accuracy  to  resolve  the  physical  behavior  of  the  weldpool  problem  at  this 
interface.  Being  too  costly  to  apply  a  global  fine  domain,  a  non-uniform  domain  is  developed 
to  utilize  finer  resolution  along  the  interface  while  still  maintaining  a  coarser  resolution  on 
the  rest  of  the  domain.  The  fast  adaptive  composite  grid  method  (FAC)  is  introduced, 
incorporating  the  concepts  of  multigrid  to  solve  the  weldpool  problem  on  this  non-uniform 
discrete  domain.  FVE  is  then  applied  to  the  conduction  equation  in  the  solid  and  the 
convection-diffusion  equation  in  the  liquid  metal  to  develop  stencil  equations  for  use  in 
FAC. 
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I. 


INTRODUCTION 


Several  practical  materiaJs  processes,  such  as  welding  and  float- zone  puriflca- 
tion,  involve  a  pool  of  molten  metal  with  a  free  surface,  with  strong  temperature  gradients 
along  the  surface.  In  certain  instances,  thermocapillary  forces  are  the  primary  drivers  of 
convection  in  the  liquid  metal.  In  cases  where  other  forces  are  stronger  overall,  the  thermo- 
capillary  forces  may  still  be  dominant  at  the  edges  of  the  molten  pool.  Strong  thermo  capillary 
convection  can  lead  to  localized  intense  heat  transfer  and  high  velocities  where  the  liquid 
free  surface  contacts  the  solid.  This  localized  heat  transfer  can  modify  the  shape  of  the 
solid-Uquid  interface  that  bounds  the  molten  pool  [Ref.  1]. 

Solving  problems  such  as  this  can  be  challenging.  One  question  that  arises  immedi¬ 
ately  is  how  to  discretize  a  continuum  problem  of  this  t5q)e.  It  is,  then  necessary  to  construct 
a  numerical  process  that  will  solve  the  resulting  discrete  system,  which  also  poses  diflRcul- 
ties.  The  purpose  of  this  work  is  twofold.  First,  the  weldpool  continuum  problem  must  be 
discretized  and  some  numerical  process  developed  to  solve  the  discrete  system.  Methods  for 
discretization  will  be  discussed,  and  then  incorporated  into  a  method  for  flnding  a  solution 
to  the  weldpool  problem.  Next,  tools  will  be  developed  that  result  from  this  discussion. 

Several  methods  for  discretizing  and  solving  the  system  will  be  discussed.  One  method 
is  the  finite  volume  element  (FVE)  method,  which  has  been  shown  to  be  an  effective  instru¬ 
ment  in  developing  discretizations,  particularly  for  problems  that  require  the  enforcement 
of  conservation  laws.  Once  discretized,  the  system  of  equations  must  be  solved.  Another 
technique,  multigrid,  has  enjoyed  remarkable  success  in  solving  certain  types  of  linear  and 
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non-lineax  partial  differential  equations  (PDEs)  in  a  uniformly  discretized  domain. 

Another  aspect  of  the  problem  results  from  examination  of  the  solid-liquid  interface. 
Near  this  boundary,  the  full  conduction-convection  solution  of  the  weldpool  problem  will 
require  extremely  fine  length  scales  to  resolve  the  physical  behavior  of  the  system.  Multigrid 
is  an  efficient  solver  on  a  grid  of  uniform  length  scale,  but  it  is  not  effective  in  cases  that 
require  fine  grid  resolution  and  accuracy  on  only  a  portion  of  the  global  domain.  The  fast 
adaptive  composite  grid  method  (FAC)  is  a  natural  extension  of  the  multigrid  methods 
applied  to  problems  discretized  on  domains  of  varying  length  scales.  Again,  FVE  serves 
admirably  as  the  tool  for  spatial  discretization. 

To  better  explain  how  the  above  techniques  can  be  applied  to  the  weldpool  prob¬ 
lem,  underlying  ideas  are  first  developed  by  applying  the  techniques  to  simple  model  prob¬ 
lems.  For  instance,  to  illustrate  the  method  of  finite  differences,  the  time  derivatives  of 
the  one-dimensional  diffusion  equation  are  discretized.  To  illustrate  the  finite  volume  ele¬ 
ment  approach,  the  potential  flow  problem  in  two-dimensional  space  is  discretized.  Having 
developed  an  understanding  of  FVE  discretization,  the  system  of  equations  governing  a  sim¬ 
plified  weldpool  problem  are  discretized,  taking  into  account  the  axisymmetric  nature  of  the 
problem. 

FAC  is  then  introduced  as  a  method  to  deal  with  the  need  for  extremely  fine  length 
scales  along  the  boundary  of  the  liquid-solid  interface.  FAC  is  an  extension  of  multigrid 
concepts.  Thus,  multigrid  is  first  developed  for  two  model  problems,  a  two-point  boundary 
value  problem  describing  the  steady-state  temperature  of  a  long  uniform  rod,  and  the  second- 
order  PDE  known  as  Poisson’s  equation.  With  an  understanding  of  multigrid  in  hand,  FAC 
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is  developed  by  building  on  concepts  from  multigrid  applied  to  a  domain  of  varying  length 
scale.  Again,  a  simplified  model  problem  is  used  to  illustrate  these  techniques.  In  this 
case,  the  model  problem  is  the  potential  fiow  problem  in  two-dimensional  space.  Finally, 
discretization  of  the  simplified  weldpool  problem  is  developed  for  a  domain  with  multiple 

length  scales. 

The  discretized  equations  developed  are  for  a  model  of  the  weldpool  problem,  while 
the  real  problem  is  far  more  complicated.  This  work  contributes  to  the  process  of  solving 
the  general  weldpool  problem  by  providing  numerical  tools  for  a  simplified  version  of  the 
problem.  Also,  complications  and  difficulties  in  the  general  problem  may  come  to  light  by 
running  the  FAC  discretization  process  on  the  simplified  problem. 
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II.  PROBLEM  STATEMENT 


A  semi-infinite  block  of  metal  is  subjected  to  a  concentrated  heat  source  ap>- 
plied  to  the  horizontal  surface  and  has  an  inviscid  nonconducting  gas  above  the  surface. 
The  result  is  a  pool  of  molten  metal  with  a  free  surface  surrounded  by  the  solid  portion 
of  the  unmelted  block  (Figure  1).  The  following  development  of  the  weldpool  problem  is  a 
recapitulation  of  Canright  [Ref.  1]. 

The  surface  tension  of  the  liquid  is  assumed  strong  enough  to  keep  the  free  surface  flat, 
but  with  surface  tension  variations  resulting  from  a  linear  dependence  on  temperature,  T. 
The  resulting  thermal  and  flow  fields  are  assumed  to  be  axisymmetric  and  steady.  The  time- 
dependent  equations  are  given  to  facilitate  a  numerical  approach  using  time-like  iterations 
to  reach  a  steady  state  solution  [Ref.  1].  The  governing  system  of  equations  is 


solid  ; 

f  = 

(11.1) 

liquid  : 

^  f = 

kV^T 

(II.2) 

—  -f-u  •  Vu  = 

-- Vp  +  i/V^u 

(11.3) 

dt 

P 

V  •  u  =  0, 

(11.4) 

with  the  conditions  at  the  boundaries  and  at  the  solid-liquid  interface  given  as 

dT 

solid  surface  z  =  0  ;  -rr  —  0 

az 

dT 

liquid  surface  z  =  0  ;  k—  =  —q{r) 

V  =  0 

du  dT 

^  dz  dr 
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This  system  is  governed  by  conservation  of  energy  in  the  solid  and  by  conservation  of  energy, 
momentum,  and  mass  in  the  pool.  In  the  system,  k  is  thermal  diffusivity,  t  is  time,  p  is 
density,  p  is  pressure,  u  is  kinematic  viscosity,  u  is  velocity  with  components  u  and  v  in 
the  r  and  2  directions  (cylindrical  coordinates),  k  is  thermal  conductivity,  q{r)  is  imposed 
surface  heat  flux  (large  at  r  =  0  and  falling  off  to  zero  at  some  small  value  of  r  such  that 
/o°°  Q{f')27rrdr  =  Q),  p,  is  viscosity,  and  7  is  the  negative  of  the  derivative  of  the  surface 
tension  with  respect  to  temperature  (7  is  assumed  to  be  constant  and  positive). 


Figure  1.  A  block  of  pure  metal  is  subjected  to  a  concentrated  surface  heating,  Q,  that  results 
in  a  molten  pool. 


The  equations  can  be  nondimensionalized.  The  main  dimensionless  parameters  are 
the  Marangoni  number  Ma  and  the  Reynolds  number  Re.  It  is  also  convenient  to  eliminate 
the  pressure  by  adopting  a  stream  function/ vorticity  formulation  for  the  flow  [Ref.  1]: 


Re  -  V  X  (u  X  (j)  j  = 

U  = 


— V  X  V  X  a; 

VxVx(|«) 

— — r  +  -— z, 
r  oz  r  or 


where  ^  is  the  axisymmetric  stream  function  and  u)  is  the  vorticity  vector  (having  only  one 
component,  in  the  9  direction).  The  flow  boundary  conditions  are 


liquid  surface  z  =  0 

:  ^  =  0 

dT 

dr 

axis  r  =  0 

:  ^  =  0 

u;  =  0 

liquid/solid  boundary  r  =  f{z,  t) 

,  d'H  ^ 

All  of  this  results  in  the  following  equations: 

solid  :  ^  =  Ma-^V^T 
at 

(II.5) 

liquid  :  ^  +  V  •  (uT)  =  Ma  ^V^T 
at 

{11.6) 

duj  „  ,  . 

—  -  V  X  (u  X  w)  = 

=  -Re  ^  V  X  V  X  cj 

(II.7) 

a,  =  VxVx  (^) 

9 

(11.8) 

where  u  =  V  x  9  =  - 

- r  H - Z 

r  az  r  ar 

(II.9) 
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with  the  conditions  at  the  boundaries  and  at  the  solid-liquid  boundary  given  as 


solid  surface  2  =  0 
liquid  surface  2  =  0 


axis  r  =  0 


dz 

dz 


=  0 

=  -q{r) 


=  0 


U 


dr 

f:=o 

or 


^  =  0 


a;  =  0. 


Refer  to  Canright  [Ref.  1]  for  a  more  in-depth  development  of  the  weldpool  problem. 

This  system  of  equations  and  systems  like  it  present  interesting  challenges  for  con¬ 
structing  numerical  solutions.  The  goal  is  to  transform  this  continuum  problem  into  a 
discrete  one  and  to  develop  a  numerical  process  to  solve  it.  The  Finite  Volume  Element 
Method  (FVE)  has  proven  to  be  a  useful  tool  in  discretizing  systems  involving  the  enforce¬ 
ment  of  conservation  laws.  Multi-level  techniques,  such  as  multigrid,  have  been  effective  in 
streamlining  the  solution  of  partial  differential  equations  (PDE’s).  The  Fast  Adaptive  Com¬ 
posite  method  (FAC)  is  designed  to  efficiently  discretize  and  find  solutions  for  PDE’s.  FAC 
methods  are  characterized  by  their  use  of  a  composite  grid,  the  union  of  various  uniform 
grids.  FAC  methods  require  discretization  of  the  PDE  on  both  the  individual  levels  and  the 
composite  grid.  FVE  is  utilized  for  this  purpose. 
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III.  DISCRETIZATION 


Numerically  solving  the  conduction  problem 


II 

{111.1) 

surface 

2  =  0  : 

K—  =  -(5(r) 

r  =  0  : 

dT  ^ 

axis 

dr 

far  boundaries  r. 

2  =  1  : 

T  =  0 

requires  discretization  of  both  the  space  and  time  derivatives,  where  the  value  of  the  unknown 
function  T  is  determined  from  a  set  of  equations  which  approximates  equation  (III.l).  In 
the  following,  finite  differences  are  used  to  discretize  the  time  derivatives,  while  the  Finite 
Volume  Element  (FVE)  method  is  used  to  discretize  the  spatial  derivatives. 

A.  FINITE  DIFFERENCES  AND  CRANK-NICOLSON 

Discretizing  in  both  time  and  space  can  be  complicated,  and  it  may  be  worthwhile  to 
develop  the  method  on  a  simple  one-dimensional  equation.  To  illustrate  the  discretization 
of  the  time  derivatives,  consider  the  one-dimensional  diffusion  equation. 


dT  d^T 
dt  dx^ 


(IIL2) 


Assuming  the  function  T  and  its  derivatives  are  single- valued,  finite  and  continuous  functions 

of  X  and  t,  then  Taylor’s  theorem  can  be  used  to  approximate  each  term  of  (IIL2).  Holding 
t  constant  (i.e.,  treating  T  as  a  function  of  x  only)  and  applying  Taylor’s  theorem  gives 
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(III,3) 


r(x  +  Ax)  =  rW  +  Ax-  +  jAx^+-Arr^  +  ... 

and 

r(x  -  Ax)  =  T(x)  -  Axf  +  i  Ax^g  -  i  Ax3g  +  . . . .  (UI.4) 

Summing  these  expressions  gives 

cP-T 

T{x  +  Ax)  +  T(x  -  Ax)  =  2T(x)  +  Ax^^  +  O(Ax^), 

at^ 

where  0(Ax'^)  denotes  the  terms  containing  fourth  and  higher  powers  of  Ax.  These  terms 
are  assumed  to  be  negligible  compared  to  the  lower  powers  of  Ax.  This  gives 

d^T  1 

^  «  25  lr(x  +  Ax)  -  2T(x)  +  r(x  -  Ax)] .  (III.5) 

Let  represent  the  temperature  at  time  nAt  and  spatial  location  kAx.  Using  this  notation, 
=  T{kAx,nAt).  Then,  in  two  dimensions,  equation  (III. 5)  is  written  as 

yr  riF„  -  22?  +  rf_, 

dx^  Ax^ 

Other  derivatives  can  be  approximated  in  similar  fashion.  For  example,  subtracting 
equation  (III.4)  from  (III.3)  and  neglecting  the  terms  of  order  Ax^  and  higher  gives 

^  5^  [r(i  +  Ax)  -  T(x  -  Ax)] .  (III.6) 

Equation  (III.  6)  is  called  a  central  difference  approximation. 

For  the  time  derivative  (holding  the  spatial  term  constant),  dT/dt  may  be  approxi¬ 
mated  in  several  ways.  The  forward  difference  is  given  as 

^  -  Tfc" 

dt  At 
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The  backward  difference  is 


^ 

dt  At 

The  central  difference  is  derived  in  a  manner  similar  to  that  leading  to  equation  (III. 6).  This 
results  in 

dT  _  2^+^  -  2]^-^ 

~  2At 

Using  the  forward  finite  difference  approximation,  the  one-dimensional  diffusion  equa¬ 
tion  (III.  2)  can  be  approximated  by 

At  Ax^ 

Letting  r  =  At/Ax^,  this  equation  can  be  written  as 

=  rT”-i  -  (1  -  +  ^^+1- 

The  unknown  can  now  be  computed  explicitly  at  the  (n  -I-  l)st  time  step  using  the 
known  values  of  the  nth  time  step. 

While  this  explicit  relation  provides  a  simple  means  to  compute  unknown  values, 
the  process  is  only  stable  for  0  <  At/(Ax)2  <  1/2  [Ref.  2].  By  imposing  a  stability  limit 
on  At/(Ax)^,  explicit  finite  difference  methods  for  the  diffusion  equation  can  give  useful 
approximations.  However,  the  requirement  that  At/(Ax)^  <  1/2  places  a  severe  restriction 
on  the  interval  size  in  the  t  direction  resulting  in  an  increase  in  computation.  A  better 
method  is  needed. 

An  implicit  formula  is  one  in  which  two  or  more  unknown  values  in  the  (n  +  l)st 
time  level  are  specified  in  terms  of  known  values  in  the  nth  time  level  by  a  single  equation. 


If  there  are  M  unknown  values  in  the  (n  +  l)st  time  level,  the  formula  must  be  applied  M 
times  across  the  length  of  the  time  level.  The  resulting  system  of  M  simultaneous  equations 
specifies  the  M  values  implicitly.  A  method  from  Crank  and  Nicolson  [Ref.  2]  combines  a 
centred  difference  approximation  to  the  time  derivative  at  the  half  time-step  (n  -f  with 
the  average  of  central  difference  approximations  to  at  times  nAt  and  (n  -1-  l)At  to 


approximate  the  equation 

Uj.  ■ 

This  equation  is  then  approximated  by 

1  (T^t;  -  m*' + g-.' , 

A!  2  I  Ai“ 


(III.9) 


(III.IO) 


This  reduces  to 


+  (2  +  27-)2^+i  -  =  rl^_j  +  (2  -  2r)T^  +  (III.ll) 


where,  again,  r  =  At/(Ax)^.  There  are  now  three  unknowns  at  the  (n  -I-  l)st  time  step 
written  in  terms  of  three  known  values  at  the  nth  time  step. 

Although  Crank  and  Nicolson’s  method  is  stable  for  all  positive  values  of  r,  large 
values  of  r,  e.g.,  40,  have  been  shown  to  introduce  oscillations  in  the  solution  [Ref.  3].  To 
combat  this  problem,  a  more  general  finite  difference  approximation  is  implemented,  using 
a  weighted  average  of  the  finite  difference  approximations  to  ^  at  the  nth  and  (n  -f-  l)st 


time  steps.  This  weighted  average  is  given  by 


_  /jTTi 


22^+1  +  yn+l 


(III.12) 


where  0  <  a  <  1.  Notice  that  a  =  0  gives  the  explicit  scheme,  a  =  1  gives  a  fully  implicit 
backward  difference  method,  and  a  =  1/2  gives  the  Crank-Nicolson  method.  This  weighted 


average  is  stable  and  convergent  for  all  r  if  ^  <  o:  <  1,  but  is  stable  and  convergent  for 
0  <  o  <  I  only  when  r  <  1/(2  —  da)  [Ref.  3]. 

Following  the  method  outlined  above,  a  general  finite  difference  approximation  for 
the  time  derivative  of  the  two-dimensional  conduction  problem  (III-l) 

^  (111.13) 

at 

is  given  by 

=  a/cV^r”^^  +  {1-  a)KV^T^.  (111.14) 

At 

In  this  approximation,  is  the  temperature  at  the  nth  time  step.  The  current  time  step 
is  designated  by  n  and  the  next  time  step  by  (n  +  1).  Combining  terms  for  the  current  and 
subsequent  time  steps  yields  the  semi-discrete  relationship. 

The  above  conditions  on  a  also  apply  here. 

B.  FINITE  VOLUME  ELEMENT  METHODS 

Having  dealt  with  the  time  discretization,  it  is  now  necessary  to  address  the  dis¬ 
cretization  of  the  spatial  derivatives.  This  development  of  discretization  for  spatial  deriva¬ 
tives  closely  follows  McCormick  [Ref.  4].  The  spatial  portion  of  the  problem  is  discretized 
using  the  Finite  Volume  Element  method  (FVE).  The  classical  finite  volume  method  (FV) 
is  commonly  used  as  a  discretization  method  for  computational  fluid  dynamics  applications 
because  of  its  ability  to  be  faithful  to  the  physics  in  general  and  conservation  in  particular. 
However,  use  of  FV  requires  a  scheme  for  approximating  certain  fluxes,  which  is  often  done  in 
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an  effective  but  rather  ad  hoc  and  restrictive  way.  Thus,  the  Finite  Volume  Element  method 
was  developed  in  an  attempt  to  use  the  flexibility  of  the  flnite  element  representation  of  the 
unknown  functions  to  create  a  more  systematic  FV  methodology  [Ref.  4] .  The  central  idea 
is  to  approximate  the  discrete  fluxes  needed  in  FV  by  replacing  the  unknown  PDE  solution 
by  a  finite  element  approximation. 


Figure  2.  Control  volume  V  (dashed  lines)  and  surface  S. 


As  an  example,  consider  the  potential  flow  problem  in  Euclidean  two-dimensional 

space: 

V  •  (pVi/))  =  rj  in  fl  (III. 16) 

Ip  =  'ipo  on  (Dirichlet) 

(pV'0)  •  n  =  ipi  on  (Neumann). 

The  domain  for  this  problem  is  the  unit  square  fi  =  [0, 1]  x  [0, 1]  G  3i^.  In  this  domain, 

and  dQs  are  the  northern  and  eastern  Dirichlet  boundaries  of  the  domain  fl,  aind 
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and  dQ.w  are  the  southern  and  western  Neumann  boundaries.  Beginning  with  a  “control 
volume”  V  in  the  domain  Vt  with  surface  S  (Figure  2),  the  potential  equation  (III.16)  is 
integrated  over  the  volume  V  to  get 


f  V  •  {pVi))dV  =  /  rjdy- 
Jv  Jv 


(III.  17) 


Using  the  Gauss  Divergence  Theorem,  the  left-hand  side  volume  integral  is  transformed  into 


Figure  3.  Volume  partition  (dashed  lines)  of  Q. 


a  surface  integral, 

/  (pVV>)  -^8=  I  T]dV.  (III.  18) 

Js  Jv 

Because  this  problem  is  two-dimensional,  “volume”  refers  to  “area”  and  “surface”  refers  to 
“boundary”.  In  other  words,  volume  integrals  in  two  dimensions  are  area  integrals,  and  the 
surface  integrals  are  line  integrals  taken  over  the  bounding  curves. 

Each  side  of  equation  (III.  18)  represents  a  flow  rate  in  mass  per  unit  time,  and 
{pVt/})  ■  h  represents  the  flux  density  across  the  surface  S.  The  integral  condition  (III. 18) 
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can  be  interpreted  as  a  conservation  law  for  the  volume  V.  Conservation  can  be  enforced 
over  the  domain  by  requiring  (III.  18)  to  hold  over  the  individual  volumes.  The  potential 
equation  (III.  16)  can  then  be  discretized  by  choosing  a  finite  set  of  volumes  that  partitions  Q 
into  volumes  (Figure  3)  with  the  integral  condition  in  (III. 18)  being  imposed  on  each  volume. 


Figure  4.  Finite  element  partition  of  the  domain 


It  is  now  necessary  to  select  the  functions  to  be  used  to  approximate  tp.  It  is  well 
known  that  piecewise  linear  functions  are  often  effective  at  approximating  functions  in  finite 
element  methods.  Thus,  ^  can  be  approximated  by  continuous,  piecewise  linear  functions. 
Then,  a  basis  —  a  set  of  functions  whose  linear  combinations  determine  the  desired  space  of 
functions  —  is  selected.  The  directions  of  interest  are  x  and  y  with  a  uniform  grid  of  step  size 
h  =  jf  applied  in  both  directions.  Each  grid  square  is  divided  into  two  triangulair  elements 
by  a  diagonal  oriented  in  the  direction  of  increasing  x  and  y  (Figure  4),  and  “hat”  functions 
are  constructed  using  these  elements.  The  FVE  method  replaces  the  exact  solution  ip*  with 
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a  finite  element  approximation  ip  expressed  as  linear  combinations  of  the  basis  functions. 
At  each  grid  point,  the  value  of  the  approximation  equals  the  weight  for  the  eissociated  hat 
function.  Thus,  ip  is  expressed  in  terms  of  a  basis: 

ij=l 

where  double  subscripts,  i,j,  varying  from  0  to  N,  are  used  to  describe  the  location  of  the 
grid  point  {xi,yj).  Then,  Uij  is  the  value  of  ip  at  the  (i,j)th  grid  point,  {xi,yj),  £ind  (pij  is 
the  “hat”  function,  a  continuous  piecewise  linear  function  associated  with  the  (i,y)th  grid 
point. 

To  simplify  notation,  the  sampled  values  of  the  unknown  ip,  currently  written  as 
an  (N  +  1)  X  {N  +  1)  two-dimensional  array,  can  be  considered  a  one-dimensional  array. 
To  accomplish  this,  the  elements  ipij  are  labeled  in  lexicographic  order.  For  example,  the 
elements  in  the  two-dimensional  case  are  written  as 

V’iV.O  V’N.l  •  •  •  'lpN,N 

V’1,0  V’1,1  •  •  •  i’l.N 

fpOfi  V’0,1  •  •  •  ‘fpO,N 

These  elements  cam  be  written  as 

(V'o.o,  V’o.l,  •  •  •  )  V’O.Af.  •  •  •,fpl,N,  •  •  •  V’AT,!)  •  •  ■  ,  • 

This  array  is  normally  considered  to  be  a  column  vector.  Using  this  notation,  ip  is  expressed 
in  terms  of  a  relabeled  basis: 

(N+ip 

ip=  (III.20) 

1=1 
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Figure  5.  Hat  function  associated  with  the  [i,j)th  grid  point. 


Substituting  (III.20)  into  (III.  18)  and  integrating  over  each  volume  I4  (following  the 
same  ordering  scheme)  for  A:  =  1, 2,  ...(N  + 1)^,  the  left  side  of  equation  (III. 18)  becomes,  for 
each  k, 

^  (N+i)" 

/  pVip  •  hdSk  =  •  nkdSk-  (III.21) 

JSk  I— I  ''Sk 

This  gives  the  system  —  /,  where 

fk=  f  vdVk  (III.22) 

JVk 

(except  on  the  boundary)  and  L  is  the  {N  +  1)  x  {N  +  1)  matrix  with  entries 

Lk,i=  f  {pkVci>i)-hkdSk.  (III.23) 

JSk 

To  correctly  implement  the  boundary  conditions,  the  system  must  be  modified.  Since 
both  Neumann  and  Dirichlet  conditions  will  be  of  interest,  both  cases  will  be  examined. 
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Neumann  conditions  are  imposed  indirectly  on  i)  by  substituting  the  flux  value  V'l  from  the 
boundary  conditions  into  the  appropriate  term  in  equation  (III. 18).  The  Neumann  conditions 
in  equation  (III.16)  give 


(III.24) 


where  Sw  and  Ss  represent  the  southern  and  western  bounding  curves  of  the  volume  where 
Neumann  conditions  are  imposed. 

By  contrast,  Dirichlet  conditions  in  equation  (III.  16)  are  imposed  directly  on  ‘tp. 
Thus,  Ip  takes  on  the  value  of  -ipo  at  each  Dirichlet  grid  point,  or  rather  each  grid  point 
on  There  are  potentially  fewer  ip  unknowns  than  equations  since  there  are  now 

volumes  that  do  not  have  unknowns.  This  is  easily  avoided  by  expanding  the  volumes  at  the 
points  neighboring  the  Dirichlet  boundary  (Figure  6). 


Figure  6.  Modified  volume  partition  of  Q.  to  accomodate  Dirichlet  points  (• ). 


In  addition,  the  intersection  of  and  d^w  and  the  intersection  of  and  d^E 
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are  treated  as  Dirichlet  grid  points.  These  points  could  also  have  been  treated  as  Neumann 
grid  points.  For  full  Neumann  problems,  the  original  partition  of  Figure  3  is  used. 

It  is  now  time  to  eveduate  the  integrals  in  (III. 22)  and  (III.23).  For  the  integral  in 
(III.23),  the  rule 

•  hdS  «  p{Pt)  V(t>  ■  ndS  (III.25) 

is  used  on  each  bounding  segment,  T,  of  the  surface.  That  is,  T  is  one  of  Sn,  Se,  Ss  or  ^ly, 
and  Pt  is  the  point  of  intersection  of  T  and  the  grid  lines  passing  through  Ny,  the  grid  point 
at  the  center  of  the  volume.  For  Neumann  boundary  segments  Ss  and  Swj  t^e  rule  is 

[  ^idS  ^  MMt)\T\,  (III-26) 

Jt 

where  Mr  is  the  midpoint  of  T  and  |r|  is  the  “surface  area,”  or  length,  of  T. 

Consider  the  discretization  on  a  uniform  mesh  with  grid  size  h  —  in  both  coor¬ 
dinates.  Returning  to  the  double  subscript  notation,  i  and  j  vary  from  0  to  —  1,  where 
0  corresponds  to  the  Neumann  boundary  grid  point,  or  node,  along  dQw  and  dQs>  and 
N  —  1  corresponds  to  the  nodes  next  to  the  Dirichlet  boundary,  and  SHe-  For  a  general 
interior  node  (z  >  0,  j  <  m  —  1),  consider  one  segment  of  the  boundary,  say  the  NW  side 
(Figure  7).  To  evaluate  the  integral 

p{Pt)  /  Vi/>  •  ndS,  (III.27) 

Jnw 

first  notice  that,  with  h  = 

•  n  =  V'l,  =  (III.28) 

since  V’  is  piecewise  linear.  Then, 
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1+1, j+1 


NW 

W  ‘  \ 

NB 

WN  1 

1  '  i 

m*  ^  *  ( 

WS  1 

SW 

^ - i 

SB 

f  VV' 

Jnw 


1-1, j-1  1#J-1 

Figure  7.  Volume  of  general  interior  grid  point  {i,j). 


V'iJ+l  -  r  (lS=-(  ^  V'ij+l  ~ 


(III.29) 


This  process  is  repeated  on  all  segments  of  the  volume  boundary  (8  Ccdculations  in  all).  The 
NW  and  NE  calculations  can  be  combined  into  a  single  North  calculation.  The  same  occurs 
for  the  South,  East  and  West  sides: 


f  Vxj}  •  ndS 

JN 

J  Vip  •  ndS 

f  Vip  •  ndS 
Je 

I  Vip-hdS  = 

Jw 


V'i+lj  ~  V’tJ 

V’i-lj  V’ij 


(111.30) 

(111.31) 

(111.32) 

(111.33) 


Then,  for  the  general  interior  grid  point,  or  node, 

p{Pt)  I  V^-nd5  = /9(Pw)^i-ij+p(F£;)V>i+ij+p(FW)V’ij+i+/’(^s)V’ij-i~5I'^»ji  (111.34) 

JSj 


where  X)  =  p{Pn)  +  p{Ps)  +  p{Pe)  +  p{Pw)- 
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The  quadrature  rule  is  used  on  the  integral  in  (III. 22).  This  results  in 

[  rjdV  ^  T]i^\Vj\,  (III.35) 

JVj 

where  rjij  is  the  value  of  r]{x,y)  associated  with  Vj  and  \Vj\  =  fy,  dV  is  the  “volume,”  or 
area,  of  Vj.  For  a  general  interior  node,  \Vj\  =  /i^.  Thus,  for  a  general  interior  node,  the 
stencil  equation  is 

^*d+2 

Pi-U  Pi+\J  tpij  =  h?v{ih,jh),  (III.36) 

Pid-^ 

where  S  signifies  the  sum  of  the  outer  stencil  entries  and  pij  =  p{ih,jh).  This  notation  for 

p{ih,jh)  will  also  be  used  in  subsequent  stencils.  In  the  stencil  notation, 

-ABC' 

DBF  'tpij, 

.G  H  1  _ 

stencil  entries  (A,  B,  ...)  indicate  the  value  to  be  applied,  with  the  position  of  the  entry 
indicating  to  which  grid  point  the  value  is  to  be  applied.  Specifically,  the  center  of  the 
stencil  is  the  current  grid  point,  with  the  surrounding  stencil  entries  corresponding  to  the 
neighboring  grid  points  in  those  directions  on  the  grid.  Thus,  the  above  stencil  has  the  value 

For  a  side  Neumann  node,  stencil  development  follows  a  similar  path.  Along  the 
boundary,  reczill  that  the  Neumann  condition  gives  'Vip-n  =■  rpi.  Integrating  on  the  Southern 
boundary  (Figure  8)  and  using  the  midpoint  rule, 

/  i>i{x,y)dS  =  hxpi{xi,0). 

J  do. 


(III.37) 


For  the  smaller  volume, 


TjdV  =  TJij\Vj\  =  T]ij  .  (III.38) 

For  a  side  Neumann  boundary  grid  point  {j  =  0  on  southern  boundary),  the  stencil  is 

5ft-i0  -E  3ft+t,o  ^,o  =  jri(ih,0)-hMih.O).  (III.39) 

For  the  corner  Neumann  node  (i  =  J  =  0),  the  stencil  is  developed  in  the  same  manner  as 


i-1,0  i,0  i+1/0 

Figure  8.  Side  Neumann  grid  point  on  the  southern  boundary  {j  =  0). 

the  side  Neumann  nodes.  The  stencil  is 

iSo.  =  ^  (lOi  (o,^)  +  V-.  .  (111.40) 

The  stencils  for  the  neighboring  Dirichlet  grid  points  (i  =  -  1  or  j  =  iV  -  1)  are 

the  same  as  at  interior  grid  points,  except  that  the  Dirichlet  points  are  eliminated  and  the 
stencil  entry  reaching  to  a  Dirichlet  node  value  is  removed  and  placed  on  the  right-hand 
side  as  a  coefficient  of  the  boundary  data.  There  are  two  possible  stencils  for  a  neighboring 
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Dirichlet  node:  the  reduced  volume  and  the  extended  volume.  For  the  neighboring  Dirichlet 
nodes  of  the  Eastern  boundary,  d^E  {i  =  N  —  <  j  <  N  —  the  reduced  volume  stencil 


h,jh)  -  (III.41) 

PN-lj-X 

where  E  is  the  sum  of  the  outer  stencil  entries  before  the  boundary  terms  are  removed.  For 
this  case,  E  =  p(l  -  ^Jh)  +  p{l  -  f,jh)  +  p{l  -  h,  {j  -  ^)h)  +  p{l  -  h,  (j  +  ^)h). 

The  stencil  for  the  same  node  in  an  extended  volume  is 


Pn-^^ 


2Pn-ij+^ 

-E 


3/j2  2 


-pNj‘^a{l,  {j  +  l)h)  -  -pNj^o{l,  {j  -  l)/i) 


where  E  is,  again,  the  sum  of  the  entries  of  the  stencil  without  the  boundary  terms  removed. 
The  extended  volume  stencil  for  the  neighboring  corner  Dirichlet  point  (i  =  N  -  1  and 
i  =  A''  —  1)  is 


Ipn-^,n-i  -  E 


i’N-ij 


2Pn-\,N-\ 


9h^  1 

-j-77(l  —  h,  i  —  h)  —  “Pat— i,Ar-i^o(l)  1  —  h) 


1  -  2/i)  -  -p^_i^^_^xpo{l  -  h,  1) 


i'0o(l  —  2/l,  1). 


C.  DISCRETIZATION  OF  THE  WELDPOOL  PROBLEM 


The  technique  used  on  the  potential  problem  (III.  16)  may  be  used  to  discretize  the 
weldpool  problem.  To  begin,  the  FVE  stencils  for  the  conduction  problem  (II.5)  are  devel- 
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Figure  9.  Toroidal  volume  for  the  conduction  problem. 

oped.  For  the  axisymmetric  conduction  equation,  a  control  volume  for  the  Ith  grid  point  is  a 
toroidal  prism.  This  volume,  VJ,  is  generated  by  taking  a  two-dimensional  volume  for  the  Ith 
grid  point  in  the  (r,  2)  plane  and  sweeping  the  volume  through  360  in  the  (p  direction  (Figure 
9).  Previously,  a  general  finite  difference  semi-discrete  approximation  for  the  equation 

^  =  Ma-^V^T 
dt 


was  found  to  be 


_  jnn 


At 


=  aMa-‘ +  (1  -  a)Ma-^V^T^. 


Combining  terms  for  the  current  and  subsequent  time  steps  yielded  the  semi-discrete  rela¬ 
tionship, 

T"*'  =  (^  +  (1  -  a)Ma-‘ T”.  (III.42) 


To  complete  the  discretization,  discrete  equations  for  the  spatial  derivatives  are  developed. 
‘  The  resulting  equations  are  assembled  into  an  algebraic  system.  Begin  by  integrating  (III.42) 


over  all  of  the  control  volumes  to  yield 


^  -  aMa-‘V2^  ^ 

Here,  V  represents  the  partition  of  the  domain  Q,  into  the  smaller  control  volumes,  Vi,  or 
V  =-{jVi.  Applying  the  Gauss  Divergence  Theorem  gives 

^  VT^+^-hdS  =  j^T^dV +{l-a)Ma-^  V'T-hdS,  (III.44) 

where  n  is  the  outward  normal. 


As  in  the  model  problem,  T  can  be  approximated  by  continuous,  piecewise  linear 
functions.  Using  cylindrical  coordinates,  a  basis  is  chosen  for  the  space  in  which  these 
functions  are  found.  Because  of  cylindrical  S5unmetry,  the  discretization  can  be  examined 
on  the  half  plane  0  =  0.  Hence,  T  is  treated  as  though  it  is  a  function  defined  on  a  two- 
dimensional  domain  in  the  (r,  z)  plane.  In  this  two-dimensional  domain,  “hat”  functions  , 
z),  are  chosen  as  the  basis  functions,  and  the  unknown  function  T  is  approximated  as 

N 

T{r,z)=  Tij<f>iAr,z),  (III.45) 

ij=0 

where  the  uniform  grid  step  size  is  h  =  1/A. 

A  change  of  notation  is  again  used.  The  sampled  values  of  the  unknown  T,  now 
written  as  an  (A -I- 1)  x  (A-F  1)  two-dimensional  array,  can  be  considered  a  one-dimensional 


array.  To  accomplish  this,  the  elements  Tij  are  labeled  in  lexicographic  order  from  /  =  1  to 
I  =  {N  1)^.  For  example,  the  elements  in  the  two-dimensional  case  are  written  as 
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These  elements  can  be  written  lexicographically  in  a  column  vector  as 


/rri  rri  rri  rrt  rri  rrt  rri  rp  rp  \i 

•  •  • )  ■‘O.AT,  1 1,0)  J  1,1)  •  •  •  •  •  ■ )  -*^^,0)  -t  yv,l)  •  •  •  i-iN,N)  ■ 


The  values  on  the  grid  Q  are  given  now  as 


Ti  =  T{ri,zi). 


Substituting  equation  (III.45)  into  equation  (III.44)  gives 


(N+iy  r  1  /■  r  1 

"  aMa-^  -hdS^  (III.46) 

(N+i)^  r  1  /•  r  1 

=  Y.  r,"  —J^(l>'^dV  +  il-a)Ma-^J^V<f)^-ndS  . 

l=X  L 

The  volume  V  of  the  domain  is  partitioned  into  I  —  {N  + 1)^  volumes  (actually,  there  are 
fewer  volumes  because  of  boundary  conditions)  where 

F  =  UF,  and  =  0  ^  ^ 

Substituting  the  partitioned  volumes  into  equation  (III.46)  results  in 

(N+lf  f  ,  (N+l)2  (AT+i)* 

E  Tr'  ^  E  /  Y  /  •  nd5,  (III.47) 

1=1  fc=l  k=l  •'^k 

(N+i)^  r  1  (N+i)^  (Ar+1)2 

=  E  J;"  ^  E  /,  +  (1  -  E  L  ■ 

1=1  fc=l  Jk=l  •'^k 


Notice  here  that  the  interior  surface  integrals  cancel  out.  The  eastern  edge  of  one  interior 
volume  coincides  with  the  western  edge  of  the  neighboring  volume.  This  edge  is  integrated 
over  twice,  but  in  opposing  directions,  resulting  in  the  cancellation. 

This  set  of  equations  can  be  written  as  an  operator  equation,  L[T^'^^]  —  /(T"),  or  as 
a  matrix  equation,  =  M2[r"].  The  operator  L  and  the  elements  of  the  matrix  Mi 
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are  given  as 


^  =  X  (s  - 

=  T:  f  't>^*'dVi-aMa-'  !  ■  ndSi,, 

At  Jvk  JSk 

while  /(T")  and  the  elements  of  M2  are  given  as 

/(T”)  =  +  (1  -  a)Ma"*V2j  T^dV 

Mfc,,  =  ^/  </>fdVk  +  {l-a)Ma-^  f  V<f>^-ndSk. 

lAZ  JVk  *'5jt 


(111.48) 

(111.49) 

(111.50) 

(111.51) 


Any  function  whose  coefficients  satisfy  the  resulting  set  of  discrete  equations  also 
satisfies  the  conservation  law  over  any  volume  made  up  of  the  union  of  control  volumes. 
However,  complications  arise  when  considering  the  boundaries.  The  Neumann  conditions 
on  the  western  and  southern  boundaries  are  incorporated  indirectly  into  T  by  substituting 
the  normal  derivative  value  into  the  appropriate  term  in  equation  (III.44).  The  Dirichlet 
conditions  on  the  northern  and  eastern  boundaries  cire  imposed  directly  on  T.  The  control 
volumes  are  extended  so  that  the  control  volumes  normedly  associated  with  the  boundary 
points  are  taken  into  account  by  the  extended  volume  of  the  neighboring  points  (Figure 
10).  The  integrals  over  the  basis  functions  are  calculated  by  representing  the  basis  functions 
as  collections  of  triangular  planes.  These  triangular  planes  are  assembled  to  form  the  hat 
function  (Figure  5).  An  interior  grid  point  has  six  triangular  planes  joined  in  this  manner, 


each  triangular  plane  corresponding  to  the  six  triangles  surrounding  the  grid  point  (Figure 
7).  Integration  results  in  the  following  stencils  for  a  general  interior  point  on  a  uniform  grid: 


32-  lie 
16 -5e 


32 -5e 
224 

32  +  5e 


16  +  56  ■ 

32  +  lie  Tij 
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Figure  10.  A  cross  section  in  the  (r,z)  plane  of  a  finite  volume  element  grid  with  the  domain 
partitioned  into  control  volumes  (dashed  lines).  Extended  control  volumes  are  used  for  grid 
points  neighboring  Dirichlet  points.  The  triangular  finite  elements  (solid  diagonal  lines)  form 
the  piecewise  linear  “hat”  function  centered  at  each  node. 


2  —  c  — 8  2  +  e 
2 


Ti 


where  e  =  ^  (r  is  the  radial  coordinate  for  the  central  point  of  the  volume).  The  stencils 
use  the  two-dimensional  labeling.  The  summations  signify  summing  over  all  “hat”  functions 
for  each  volume.  Also,  the  volume  and  surface  integrals  are  in  cylindrical  coordinates.  This 


means  that  the  integrals  are  as  follows: 


f  dV  =  f  rdrdOdz  =  2ir  f  rdrdz 
Jv  Jv  Ja 

j  dS  =  j  rdOdl  =  27r  ^  rdl. 

In  each  case,  the  27r  resulting  from  the  integration  in  the  9  direction  has  been  factored  out. 
The  discretization  of  the  convection-diffusion  equation, 

^  +  V  •  (uT)  =  Ma-^V^T, 
at 
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has  the  same  form  as  the  conduction  equation  for  the  solid  except  for  the  term 


V  •  {uT)dV. 

Using  the  Gauss  Divergence  Theorem,  this  portion  of  the  equation  is  also  discretized.  By 
cylindrical  symmetry, 


Substituting  the  piecewise  linear  “hat”  functions  for  the  unknowns  and  T  into  the  above 
integral  results  in  discrete  equations.  Evaluating  and  summing  the  integrals  and  using  stencil 
notation  gives 


30 


Discretization  of  the  vorticity  (equation  II.7)  follows  the  same  scheme.  Integrating 


(III.52) 


the  vorticity  equation  over  a  control  volume  results  in 

^  j  Cj  •  6dV  —  J  V  X  (u  X  (l))dV  =  —Re~^  ^  CjdV.  (III. 52) 

The  integration  over  the  volume  equates  to  an  integral  over  the  cross  sectional  area  with  Q 
acting  as  the  normal  vector.  Applying  Stokes’  theorem  to  (III.52)  results  in 

J  J  (j^drdz  —  j>  t  •  (u  X  (jj)dl  =  —Re~^  ji  i  ■  (V  x  co)dl.  (III. 53) 

Again,  a  uniform  grid  with  step  size  h  in  both  the  r  and  z  directions  is  applied.  Each  square 
of  the  grid  is  divided  into  two  triangular  elements  by  a  diagonal  in  the  direction  of  increasing 
r  and  z  (Figure  10)  and  linear  interpolation  is  used  over  each  triangular  element.  The  control 
volume  cross  sections  are  square  with  side  length  h  and  centered  on  a  grid  point,  except  for 
the  control  volumes  surrounding  side  and  corner  Neumann  grid  points  and  side  and  corner 
Dirichlet  boundary  points  (again,  see  Figure  10). 

The  area  integrals  are  taken  over  six  triangular  regions,  while  the  line  integrals  are 
taken  over  four  line  segments,  each  line  segment  containing  halves  in  two  different  elements. 


Recall  that 


1  ,  Id'^ , 

— ^r  +  — ^z, 
r  oz  r  or 


(III.54) 


and  note  that  the  vorticity  vector,  u,  has  only  one  component,  the  component  in  the  9 
direction.  Using  the  expression  for  u  from  above  results  in 


,  u}  d'if  ^  u)  d'^ , 

u  X  a;  =  — — r - — z. 

r  or  r  oz 


The  integrals  in  (III.53)  then  reduce  to 

d  r  r  r  u)  ,  f  d'i!  u  ,  r  d'i!  u  ,  r  uj 

—  /  /  uidrdz  +  /  — - dr  —  /  — dr  +  /  — - dz  —  /  — - 

dt  J  Ja  Jn  dr  r  Js  dr  r  Jw  dz  r  Je  dz  r 


=  Re 


rpr-f^-^^dz+f?P-^-dz). 

Js  oz  Jw  or  r  Je  or  r  I 


After  substituting  the  piecewise  linear  element  representation  of  the  unknowns  into 
the  above  integrals,  the  discrete  equations  for  an  interior  grid  point  on  a  uniform  grid  are 
given  in  stencil  form  as 


I  J'l 

— — —  2  14  2  u 

dt  24  1  1  2  ) 


2:^1+ 

-1-6+ 


-2  l-ie+ 


-1+^6+ 


^  H-C+  -e+  +  e~  l-e“  Mk 

V  -1-K  J 


1  +  fe-  -2  ^ 


-1 


1-  ^ 


— 1  +  e" 
2  *  1  — 


^  u 


=  r(ife-‘)  l-fe+  -4 -§(/+-£-)  1  +  |£- 


1  +  K 


where  c  =  h/r,  e'^  =  e/(l  —  e/2),  e  =  e/(l  +  e/2),  1"^  =  1/(1  —  e/2),  and  1  =  1/(1  +  e/2). 
The  last  equation  to  be  considered  for  discretization  is  the  stream  function  (II.8), 


o)  =  V  X  u, 


(111.55) 


where 


u  =  Vx(^-j0. 


Integrating  over  the  control  volumes  and  applying  Stokes  Theorem  to  the  right  side  of  equa¬ 


tion  (III. 55)  gives 


j  j  Ljdrdz  —  ^ 


(III.56) 
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where  i  is,  again,  the  unit  vector  tangent  to  the  curve  C.  With  u  as  above  (equation  III.  54), 
these  integrals  reduce  to 


uidrdz 


-L 


- — dr-\-  /  - — dT  + 
N  oz  r  Js  dz  r 


f  1 
Is  dz  r 


[ 

Iw  dr  r 


I 

Ie  dr  r 


(111.57) 


to  resolve  the  weldpool  problem  more  efficiently. 
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IV.  FAST  ADAPTIVE  COMPOSITE  METHODS 


This  chapter  will  introduce  the  concepts  of  fast  adaptive  composite  (FAC) 
grid  methods.  FAC  methods  are  used  to  discretize  and  find  solutions  for  partial  differential 
equations  (PDFs)  on  a  composite  grid.  A  composite  grid  is  formed  by  the  union  of  a  nested 
sequence  of  uniform  grids  of  varying  size  and  scale.  A  patch,  or  rectangular  uniform  grid, 
can  be  aligned  with  coarser  grids  and  can  be  specified  by  its  relative  location  on  the  coarser 
grid,  its  dimensions,  and  its  mesh  size  (Figure  11).  The  PDE  is  discretized  and  solved  on 
this  composite  grid.  However,  the  actual  computations  take  place  on  the  uniform  subgrids. 
In  this  manner,  FAC  maintains  the  advantages  of  uniform  grid  discretizations  while  allowing 
effective  adaptation  to  the  local  phenomena.  FAC  methods  require  discretization  of  the  PDE 
on  both  the  individual  levels  and  the  composite  grid.  The  approach  used  here,  the  finite 
volume  element  (FVE)  method,  combines  the  finite  volume  and  finite  element  methods  [Ref. 

4]. 

Many  physical  models  encountered  provide  excellent  opportunities  to  employ 
FAC.  Often,  physical  models  require  higher  accuracy  only  in  limited  regions  of  a  domain. 
The  physical  model  could  be  resolved  on  a  global  uniform  grid  with  the  smallest  required 
mesh  size,  but  computational  power  is  rarely  unlimited,  and  some  method  for  locally  adapt¬ 
ing  resolution  and  accuracy  is  needed.  It  can  also  be  an  advantage  to  introduce  as  little 
nonuniformity  into  the  method  as  possible.  This  is  accomplished  using  composite  grids.  For 
an  example  of  a  coarse  grid,  consider  the  two-level  grid  in  Figure  12.  The  coarse  grid  is 
defined  on  a  domain  where  2h  indicates  the  grid  spacing.  There  are  several  distinct 
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Figure  11.  A  composite  grid  with  its  fine  grid  patch. 


types  of  gridpoints,  such  as  regular  grid,  coarse  grid  interface,  fine  grid  interface,  and  slave 
points.  The  fine  grid  points  can  be  treated  as  regular  interior  points  of  the  refinement  patch, 
while  coarse  grid  interface  points  and  slave  interface  points  can  be  taken  as  boundary  points 
for  the  fine  grid  patch.  Since  the  functions  on  all  of  the  grids  are  piecewise  linear,  the  slave 
interface  points  can  be  calculated  as  a  linear  interpolation  of  neighboring  coarse  grid  in¬ 
terface  points.  Multilevel  methods  can  then  be  implemented.  These  methods  make  use  of 
fully  overlapping  grids  of  different  scales  and  correction  of  the  coarse  grid  equations  by  the 
composite  grid  residuals  to  approximate  the  PDE  solution.  This  is  accomplished  by  solving 
the  system  on  the  fine  grid  level  and  then  correcting  the  coarse  grid  equations  with  the  fine 
grid  approximations. 
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^  fin*  grid  int*rf*o* 


*l«v*  point*  conr**  grid  lnt«rf*o* 


Figure  12.  Composite  Grid  Interface  Points. 

A.  FUNDAMENTALS  OF  MULTIGRID 

Before  entering  into  a  development  of  FAC,  the  concepts  of  multigrid  should  be  in¬ 
troduced.  Multigrid  is  a  basic  component  of  FAC  because  of  its  apparent  effectiveness  as 
an  iterative  method.  It  is  of  interest  to  determine  how  multigrid  behaves  as  a  solver  of 
discretizations  on  uniform  grids,  because  multigrid  will  be  used  for  composite  grid  equations 
restricted  to  the  global  and  local  subgrids  [Ref.  4]. 

Multigrid  methods  were  first  developed  to  solve  boundary  value  problems  posed  on 
spatial  domains.  These  boundary  value  problems  can  be  made  discrete  by  selecting  a  set  of 
grid  points  in  a  domain.  The  discrete  problem  that  results  is  a  system  of  algebraic  equations 
associated  with  the  chosen  grid  points.  The  simplicity  of  these  boundary  value  problems 
provides  a  good  introduction  to  multigrid. 


1.  Two  Model  Problems 


Consider  the  two-point  boundary  value  problem  describing  the  steady-state  temper¬ 
ature  distribution  in  a  long  uniform  rod.  The  second-order  ordinary  differential  equation 
is 

—u"{x)  +  au{x)  =  /(x),  0  <  X  <  1,  cr  >  0,  (IV.  1) 

where  u(0)  =  u(l)  =  0.  The  simplest  method  for  solving  this  equation  numerically  is  the 
finite  difference  method.  The  domain  of  the  problem  is  partitioned  into  N  subintervals  in  a 
step  size  oih  =  1/N.  A  grid  is  established  with  points  Xj  =  jh  (Figure  13). 
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Figure  13.  One- dimensional  grid  on  the  interval  0  <  x  <  1  vnth  grid  spacing  h  =  1/N  and 
grid  point  Xj  =  jh  for  0  <  j  <  N . 


At  each  of  the  N  —  1  interior  grid  points,  the  differential  equation  (IV.l)  is  replaced  by 
a  finite  difference  approximation.  Since  the  exact  solution  is  unknown,  the  approximation  Vj 
is  introduced.  The  approximate  solution  is  represented  by  the  vector  v  =  (vi,V2, ,  uat-i), 
whose  components  satisfy  the  TV  —  1  linear  equations  of  the  finite  difference  approximation. 


-Vj+i  +  2vj  -  Vj-i 

h? 


+  avj  =  /(xj), 


1  <  i  <  -V  - 1, 


(IV.2) 
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with  the  boundary  conditions  vq  =  vjv  =  0-  Representing  this  system  of  linear  equations  in 
matrix  form  gives 

"  2  +  ah?  —1  /i 

— X  2,  ah?  — 1  ' 

i_  ...  •  ^  3^ 

h2  ...  .  • 

.  .  -1 

-1  2  +  ahi?  J  [  vn-i  J  L  /n-1  . 

or  simply  Av  =  f.  The  matrix  A  is  tridiagonal,  symmetric  positive  definite,  sparse  (the 

matrix  elements  are  predominantly  zero),  and  has  dimension  {N  —  1)  x  {N  —  1). 

For  the  two-dimensional  case,  consider  Poisson’s  equation,  a  second-order  partial 
differential  equation  given  by 

-Uxx  -  Uyy  =  f{x,  y),  0  <  X  <  1,  0  <  y  <  1,  (IV.4) 

with  u  =  0  along  the  boundary  of  the  unit  square  (Figure  14). 


Figure  14.  Two-dimensional  grid  on  the  unit  square. 


As  before,  Poisson’s  equation  can  be  cast  in  a  discrete  form  using  the  grid  points 
{xi,yj)  =  {ih,jh),  where  h  =  1/N  is  the  width  of  the  subintervals.  Replacing  the  derivatives 
of  Poisson’s  equation  (IV.4)  with  second-order  finite  differences  and  the  exact  solution  with 
the  approximation  Vij  gives 


+  2Vij  -  Vj-ij  -Vjj+i  +  2Vij  VjJ-l  _  s 

fi2 


(IV.5) 


with  Vij  =  0  along  the  boundary  [i  =  0  or  i  =  N  oi  j  =  0  ot  j  =  N)  and  1  <  i  <  A  —  1  and 
I  <  j  <  N  -  1.  There  are  now  {N  -  1)^  interior  grid  points  and  unknowns.  By  arranging 
the  lines  of  constant  i  in  lexicographic  order,  the  unknowns  of  the  ith  row  of  the  grid  can 
be  collected  in  the  vector  vj  =  (I’i,!,  Ui,2>  •  •  •  i  Vi^N-i)'^  for  1  <  i  <  A  —  1.  When  f{xi,  yj)  is 
ordered  in  the  same  manner,  the  system  of  equations  (IV.5)  can  then  be  written  in  matrix 


form  as 


(IV.6) 


This  system  is  symmetric,  block  tridiagonal,  and  sparse.  Each  block  on  the  diagonal 
is  an  (A  -  1)  X  (A  -  1)  tridiagonal  matrix,  while  each  off-diagonal  block  is  the  negative  of 
the  (A  —  1)  X  (A  —  1)  identity  matrix  I. 

2.  The  Residual  Equation 

Remember  that  a  system  of  linear  equations,  such  as  equations  (IV.2)  and  (IV.5), 
can  be  written  as 

Au  =  f,  (IV.7) 
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with  u  being  the  exact  solution.  Let  v  represent  an  approximation  to  the  exact  solution. 
The  exact  solution  is  unknown,  while  the  approximation  v,  which  is  often  computed  using 
an  iterative  method,  is  known.  One  important  measure  of  v  as  an  approximation  of  u  is 
algebraic  error,  given  by 

e  =  u  -  V.  (IV.8) 

However,  since  the  exact  solution  is  still  not  known,  the  error  is  also  unknown.  Still,  one 
measure  of  how  well  v  approximates  u  is  the  residual,  given  by 

r^f-Av.  (IV.9) 

The  residual  is  the  amount  by  which  the  approximation  fails  to  satisfy  the  original  problem 
in  (IV. 7).  Rewriting  the  residual,  (IV.9),  as 

^v  =  f-r  (IV.IO) 

and  subtracting  this  new  version  from  the  original  system,  equation  (IV.7),  produces  the 
residual  equation 

=  r.  (IV.ll) 

Assuming  that  the  approximation  v  has  been  computed  by  some  method,  the  residual  r  can 
be  easily  computed  using  r  =  f  —  Av.  To  improve  the  approximation,  solve  the  residual 
equation  (IV.ll)  for  e  and  compute  a  new  approximation  using  the  definition  of  the  error 
(IV.8).  This  is  the  idea  of  residual  correction  and  will  prove  to  be  very  important  [Ref.  5]. 

The  residual  equation  shows  that  the  error  satisfies  the  same  set  of  equations  as  the 
desired  exact  solution  when  f  is  replaced  by  r.  This  suggests  that  performing  some  iterative 
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method  on  the  original  equation  =  f  with  an  arbitrary  guess  v  is  equivalent  to  performing 
the  same  iterative  method  on  the  residual  equation  with  the  specific  initial  guess  e  =  0. 

3.  Solving  Systems  of  Linear  Equations 

There  are  several  existing  methods  for  solving  sparse  systems  of  linear  equations 
like  equations  (IV.2)  and  (IV.4).  Direct  methods  solve  the  system  of  equations  in  a  known 
number  of  arithmetic  operations.  Several  types  of  errors  can  result  from  using  direct  methods. 
Algebraic  errors  in  the  solution  result  from  rounding  errors  introduced  by  the  computation. 
Discretization  error  is  the  difference  between  the  exact  solution  of  the  PDE  and  the  exact 
solution  of  the  difference  equations  being  used  to  approximate  the  PDE.  The  magnitude 
of  the  discretization  error  at  a  grid  point  depends  on  the  distances  between  consecutive 
discrete  grid  points  and  on  the  number  of  terms  in  the  truncated  series  of  differences  used 
to  approximate  the  derivatives. 

Direct  methods  tend  to  be  elimination  methods,  such  as  Gaussian  elimination,  or 
factorization  methods,  such  as  LU  decomposition,  QR  factorization,  or  Singular  Value  De¬ 
composition  (SVD).  As  seen  from  the  systems  above,  most  difference  approximations  result 
in  sparse  matrices,  with  the  number  of  zero  entries  greatly  outnumbering  the  non- zero  entries. 
For  this  type  of  matrix,  a  method  like  Gaussian  elimination  proves  inefficient,  because  the 
procedure  fills  in  the  zero  elements  with  non-zero  elements  during  the  elimination,  thereby 
increasing  the  number  of  non-zero  elements  that  need  to  be  stored  in  the  computer  and 
used  in  later  eliminations.  As  a  result,  direct  methods  can  be  quick  and  accurate  for  certain 
matrices  with  special  properties  but  prove  inefficient  with  elliptic  equations  [Ref.  3]. 
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Iterative,  or  relaxation,  methods  are  such  that  an  initial  approximation  is  used  to 
calculate  a  second  approximation  which  is  used  to  calculate  a  third,  and  so  on.  Iterative 
methods  are  said  to  be  convergent  when  the  differences  between  successive  approximations 
tends  to  zero  as  the  number  of  iterations  increases.  Iterative  methods  have  been  much 
employed  to  solve  vlu  =  f  because  they  take  advantage  of  the  sparcity  of  A.  These  relaxation 
methods  are  easy  to  implement  and  can  be  successful  in  determining  an  exact  solution  for 
more  general  linear  systems  than  the  direct  methods.  However,  experiments  have  shown  that 
these  methods  generally  work  very  well  only  for  the  first  several  iterations.  Convergence  then 
slows  8ind  the  method  appears  to  stall. 

Iterative  methods,  such  as  the  Jacobi  and  Gauss-Seidel  methods,  begin  with  an  initial 
guess.  Valuable  insight  can  be  gained  by  using  an  initial  guess  consisting  of  Fourier  modes: 

j,  0<j<iV,  1<A:<A^-1.  (IV.12) 

The  integer  k  is  the  wavenumber  and  indicates  the  number  of  half  sine  waves  which  make 
up  the  vector  v  on  the  domain.  Small  values  of  k  correspond  to  long,  smooth  waves.  Large 
values  correspond  to  oscillatory  waves.  The  rapid  decrease  in  error  during  the  first  several 
applications  of  an  iterative  method  results  from  efficient  elimination  of  the  oscillatory  modes 
of  this  error.  Unfortunately,  the  iteration  is  much  less  effective  in  reducing  any  remaining 
smooth  components.  Most  relaxation  methods  possess  this  property,  known  as  the  smoothing 
property  [Ref.  5].  The  issue  is  whether  basic  iterative  methods  can  be  adapted  to  make  them 
effective  on  all  error  components. 

One  way  to  improve  an  iterative  method  is  to  start  with  a  good  initial  guess.  By 
performing  some  preliminary  iterations  on  a  coarse  grid,  an  improved  initial  guess  can  be 


Vj  =  sin 


jkiT 

N 
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derived  and  used  on  the  original  fine  grid.  This  coarse  grid  is  a  grid  (Figure  14)  with  a  larger 
interval  between  grid  points.  Performing  these  iterations,  or  relaxing,  on  a  coarse  grid  is  less 
expensive  since  there  are  fewer  unknowns  to  be  updated. 

With  this  is  mind,  recall  that  most  relaxation  schemes  are  effective  at  eliminating 
the  oscillatory  components  of  error  but  prove  futile  against  the  smooth  components.  Now, 
consider  these  remaining  smooth  components  on  a  coarser  grid.  A  smooth  curve  on  a  fine 
grid  can  appear  more  oscillatory  on  a  coarse  grid.  Note  in  Figure  15  that  the  grid  points 
on  the  coarse  grid  (bottom  grid)  correspond  to  the  even-numbered  grid  points  of  the  fine 
grid  (top  grid).  Consider  the  A:th  mode  of  the  fine  grid  evaluated  at  the  even-numbered  grid 
points.  If  1  <  /:  <  N/2,  then  the  components  of  the  A:th  mode  on  the  fine  grid  can  be  written 
as 


'2jkTT 
.  N 


sin 


'  jk-K ' 


=  w 


2h 

kj' 


(IV.  13) 


The  superscripts  on  w  indicate  on  which  grid  the  vector  is  defined  by  giving  the  grid 
spacing  on  the  grid.  The  vector  is  then  defined  on  the  fine  grid,  while  the  vector  is 
defined  on  the  coarse  grid. 

Thus,  the  A:th  mode  on  the  fine  grid  becomes  the  kth  mode  on  the  coarse  grid.  The 
A:th  mode  also  becomes  more  oscillatory  in  passing  from  the  fine  grid  to  the  coarse  grid. 
FVom  figure  15,  A:  =  4  on  the  grid  where  W  =  12isl/3of  the  highest  possible  frequency, 
while  A;  =  4  on  the  grid  where  W  =  6  is  2/3  of  the  highest  possible,  making  the  wave  more 
oscillatory  on  the  coarse  grid.  This  is  true  only  for  1  <  A:  <  N/2.  If  k  =  N/2,  the  A:th  mode 
on  the  fine  grid  is  the  zero  vector  on  the  coarse  grid.  And  for  k  >  N/2,  the  kth  mode  on  the 
fine  grid  becomes  the  (N  —  A:)th  mode  on  the  coarse  grid.  When  this  phenomenon,  known 
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k  s  4  wave  on  N  s  6  grid 


Figure  15.  A  wave  with  wdvcuuTtibev  k  =  4  oti  a  fine  grid  (^N  —  12)  appears  vnore  oscillatory 
on  a  coarse  grid  {N  =  6). 

as  aliasing  ,  occurs,  oscillatory  modes  on  the  fine  grid  are  misrepresented  as  smooth  modes 
on  the  coarse  grid. 

All  of  this  suggests  that,  once  convergence  slows  on  the  fine  grid,  the  residual  may 
be  computed  and  the  residual  equation  can  be  relaxed  on  a  coarser  grid  to  obtain  an  ap¬ 
proximation  of  the  error.  This  error  can  then  be  returned  to  the  fine  grid  to  correct  the 
original  approximation.  This  method  for  improving  iterative  approximations,  called  coarse 
grid  correction,  leaves  some  issues  unresolved,  such  as  how  to  transfer  vectors  from  a  coarse 
grid  to  a  fine  grid  (and  back  again),  how  to  compute  the  residual  on  the  fine  grid,  and  how 
to  compute  an  error  estimate  on  the  coarse  grid. 

4.  Inter  grid  Transfers 

Before  an  iterative  method  can  be  used  on  coarser  grids,  a  mechanism  is  needed  for 
tranferring  vectors  between  grids.  A  discussion  of  intergrid  transfers  will  deal  only  with  the 
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case  in  which  the  coarse  grid  has  twice  the  grid  spacing  of  the  next  finest  grid,  since  there 
seems  to  be  no  advantage  in  using  grid  spacing  with  ratios  other  than  two. 

Transferring  a  vector  from  a  coarse  grid  to  a  fine  grid  is  called  interpolation,  or 
prolongation.  Although  many  methods  exist,  for  most  multigrid  purposes  linear  interpolation 
is  quite  effective.  The  linear  interpolation  operator  will  be  denoted  as  Z^/,.  The  operator 
takes  coarse  grid  vectors  and  produces  fine  grid  vectors.  With  as  the  coarse  grid  vector 
and  v^  as  the  fine  grid  vector,  the  notation  transfers  the  vectors  such  that 


— 


,2h 


.  N 

}  0<j<j-l. 


4+1  =  J 

This  transfer  notation  shows  that,  at  even-numbered  fine  grid  points,  the  values  of  the  vector 

are  transferred  directly  from  the  coarse  to  the  fine  grid.  At  odd-numbered  fine  grid  points, 

the  value  of  the  grid  point  on  the  fine  grid  is  the  average  of  the  adjacent  coarse  grid  values. 

For  a  two-dimensional  problem,  the  interpolation  operator  is  defined  in  a  similar 

manner.  According  to  the  notation  =  v'*,  the  components  of  are  given  by 

,2/1 

Oh\ 

N 

}  0<i,j<j-l. 


— 

4+1, 2j 

^2i,2j+l 

’2i+l,2j+l 

—  X  ( -4-  -L 

Transferring  a  vector  from  a  fine  to  a  comse  grid  is  called  restriction.  One  restriction 
operator,  injection,  is  defined  by  =  v^^,  where  the  components  of  v^^  are  given  by 


„,2h  _  h 
Vj  -V^j. 
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With  injection,  the  coarse  grid  vector  takes  its  value  directly  from  the  corresponding  fine 
grid  point. 

Another  restriction  operator,  known  as  full  weighting,  is  defined  by 
where  the  values  of  the  coarse  grid  vector  are  a  weighted  average  of  values  at  neighboring 
fine  grid  points,  or 

The  two-dimensional  full  weighting  operator  is  defined  similarly,  where  the  values  of  the 
coarse  grid  vector  take  on  a  weighted  average  of  values  at  neighboring  fine  grid  points. 
Thus, 

^id  =  W  [’^2i-l,2i-l  +  V2i-l,2j+l  +  ^2i+l,2j-l  +  '^2i+l,2j+l 
+2  (^2i,2i-l  '^2»,2ji+l  ^2»-l,2j  +  ^2i+l,2j) 

+'^^2i,2i]  >  1  <  ^  y  ~  1* 

5.  Coarse  Grid  Correction 

A  coarse  grid  correction  scheme  may  now  be  introduced.  Coarse-grid  correction  is 
described  as  a  two-grid  process  as  follows: 

•Relax  i/i  times  on  A'*u''  =  f  on  O'*  with  an  initial  approximation  v'*. 

•Compute  the  residual  on  the  fine  grid,  r'*  =  f -  Av'*,  transfer  the  residual  from  the 
fine  grid  to  the  coarse  grid,  and  solve  the  residual  equation,  =  r^'*, 

on  the  coarse  grid  with  an  initial  guess  =  0. 

•Correct  the  approximation  obtained  on  (v'*  in  the  first  step)  with  the  error  estimate 
obtained  on  v'*  •<-  v'*  -I- 

•Relax  1/2  times  on  A'‘u'*  =  f  on  Q'*  with  the  improved  initial  approximation 
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The  integers  Vi  and  i>2  are  parameters  in  the  correction  scheme  which  control  the  number  of 
relaxation  sweeps  before  and  after  the  coarse  grid  correction.  Experiments  have  shown  that 
effective  choices  for  Ui  and  1^2  are  often  1,  2,  or  3. 

is  the  coarse  grid  representation  of  the  original  matrix  A.  The  original  problem 
could  be  directly  discretized  on  using  the  same  techniques  used  to  generate  A^,  the  fine 
grid  operator,  but  at  a  spacing  of  2h,  resulting  in  A^^.  An  alternate  approach  is  to  note  that 
if  the  error  on  the  fine  grid,  e*  =  —  v^,  lies  entirely  within  the  range  of  interpolation, 

denoted  by  R(l2h)^  then  the  error  can  be  represented  by  for  some  coarse  grid 

vector  The  residual  equation  on  can  then  be  written  as 

Applying  the  restriction  operator  to  the  equation 

=  r* 

gives 

T2h  Ah  rh  --2/i  _  r2h^h 

4  A  Yj/iU  -  4  r 
which  leads  to  the  natural  definition 

The  coarse  grid  equation  becomes 

j^2k^2h  ^  j.2h 

This  argument  is  based  on  the  assumption  that  the  error  e*  is  entirely  within  the  range  of 
interpolation.  This  may  not  always  be  the  case.  If  the  error  were  entirely  within  this  range 
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of  interpolation,  then  solving  the  residual  equation  on  the  coarse  grid  would  give  an  exact 
solution.  Regardless,  this  approach  gives  a  definition  for  the  coarse  grid  operator  This 
definition  is  known  as  the  Galerkin  condition.  This  coarse  grid  operator  may  not  give  the 
same  as  discretizing  at  twice  the  mesh  spacing.  However,  the  Galerkin  condition  may 
be  used  to  prove  such  things  as  convergence  theorems,  making  this  condition  an  effective 
theoretical  tool. 

In  the  above  coarse  grid  scheme,  there  is  an  efficient  way  of  solving  the  coarse  grid 
problem  The  coarse  grid  problem  is  not  much  different  from  the  original 

problem,  and  can  be  solved  by  the  same  coarse  grid  correction  scheme.  Applying  the  coarse 
grid  correction  scheme  to  the  residual  equation  on  the  coarse  grid  requires  a  move  to  the 
next  coarsest  grid,  for  the  correction  step.  This  process  is  repeated  on  successively 
coarser  grids  until  a  direct  solution  to  the  residual  equation  is  economical.  The  resulting 
algorithm  is  applied  recursively  down  to  the  coarsest  grid,  sometimes  a  single  interior  grid 
point.  Corrections  are  then  applied  successively  through  finer  and  finer  grids,  returning 
eventually  to  the  finest  grid.  This  algorithm  is  referred  to  as  a  V-cycle  ,  named  for  the 
pictoral  representation  of  the  schedule  of  grid  visits  (see  Figure  16).  It  is  defined  by  ■<— 
aJid  is  summarized  here  using  a  compact  recursive  definition: 

•Relax  ui  times  on  A'*u'*  =  f  on  n'‘  with  an  initial  approximation  v'*. 

•If  fl'*  is  the  coarsest  grid,  then  proceed  to  the  last  step.  Otherwise,  compute  the 
residual  on  the  fine  grid,  r'*  =  f'*  -  Av'*,  transfer  the  residual  from  the  fine  grid  to 
the  next  coarsest  grid,  r^'*  =  and  set  the  initial  guess  for  to  zero,  ^  0. 

Return  to  the  first  step,  r^'*). 

•Correct  v'*  -f-  v'"  + 

•Relax  1/2  times  on  the  original  equation,  =  f'*,  with  an  initial  guess  v*. 
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Recall  that  one  way  to  improve  iterative  methods  is  to  start  with  a  good  initial  guess. 
The  V-cycle  uses  only  coarse  grid  correction  in  solving  the  residual  equation  and  does  not 
make  use  of  any  method  to  improve  the  initial  guess.  The  method  of  nested  iteration  uses 
coarse  grids  to  obtain  improved  initial  guesses  for  fine  grid  problems.  In  order  to  obtain 
an  improved  initial  guess  for  the  first  fine  grid  relaxation  of  the  V-cycle,  nested  iteration 
suggests  solving  this  problem  on  the  next  coarsest  grid.  To  obtain  an  improved  initial  guess 
on  this  coarse  grid,  say  it  is  necessary  to  solve  the  problem  on  O'**.  This  pattern  results 
in  a  recursive  path  that  leads  to  the  coarsest  grid. 

The  full  multigrid  V-cycle  (FMV)  joins  nested  iteration  with  the  V-cycle,  resulting  in 
a  simple  but  powerful  algorithm.  Each  V-cycle  is  preceded  by  a  smaller  V-cycle  that  provides 
the  best  initial  guess  possible.  The  full  multigrid  V-cycle,  defined  by  v^  <— 
is  summarized  here; 
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•Initialize  f'*,  r^'*,  r'^,  . . v'*,  ...  to  zero. 

•Relax  on  the  residual  equation  on  the  coarsest  grid  (or  solve  directly  if  the  coarsest 
grid  is  small) . 

•Correct  and  perform  a  V-cycle  using  the  grid  as  the  finest 

grid  and  fl®*  as  the  coarsest  grid.  This  is  abbreviated  as  •<—  MV’^^(e^,r^). 

•Correct  using  the  result  from  the  previous  V-cycle,  -1- 

•v^  •<—  -I-  /2/,V^^. 

•v'‘ ^  MV'‘(v^f'*). 

The  last  step  in  the  above  procedure  is  the  original  V-cycle  discussed  earlier,  but  is  preceded 
by  several  smaller  V-cycles  (Figure  16)  in  order  to  make  an  improved  initial  guess. 

The  full  multigrid  V-cycle  brings  together  many  ideas  and  techniques,  that,  when 
taken  alone,  can  have  serious  limitations.  Full  multigrid  integrates  these  ideas  so  that  they 
may  work  together  to  remove  these  limitations. 

B.  THE  FAST  ADAPTIVE  COMPOSITE  GRID  METHOD 

The  fast  adaptive  composite  grid  method  (FAC)  is  a  natural  extension  of  conven¬ 
tional  multigrid  methods  applied  to  problems  discretized  on  composite  grids.  Multigrid  was 
developed  as  a  global  grid  solver  but  is  not  efficient  for  cases  where  fine  grid  resolution  and 
accuracy  are  needed  only  in  a  local  subdomain.  While  a  certain  amount  of  accuracy  may  be 
suitable  to  correctly  solve  the  physical  system  in  most  of  the  global  domain,  greater  accu¬ 
racy  is  most  often  needed  in  a  given  local  region.  Obtaining  the  desired  accuracy  may  take 
a  resolution  with  a  mesh  size  h,  which  is  smaller  than  that  in  the  rest  of  the  global  domain. 
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If  a  global  fine  grid  were  incorporated,  unnecessary  computations  would  result  outside  the 
local  region,  often  at  a  great  computational  cost  [Ref.  4]. 


Figure  17.  A  global  domain  0.  Q2  is  the  local  domain  with  a  smaller  grid  spacing  while  f2i 
is  that  portion  of  the  global  domain  ft  that  retains  the  original  coarse  grid  spacing. 

To  arrive  at  FAC,  the  desired  method,  ideas  will  be  drawn  from  several  simpler,  less 
efficient  methods.  It  should  be  noted  here  that  this  development  will  follow  closely  with  that 
of  reference  [Ref.  4]. 

1.  Local  Multigrid 

A  natural  way  to  avoid  using  a  global  fine  grid  without  sacrificing  the  necessary  fine 
grid  resolution  in  the  local  subdomain  is  to  modify  the  standard  multigrid  scheme  by  simply 
suppressing  relaxation  outside  of  the  subdomain.  Let  represent  a  restricted  relaxation 
process  on  the  local  subdomain  SI2  (Figure  17).  Let  the  decomposition  of  be  represented 
by 
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where  U2  represents  the  unknown  values  of  u  found  on  the  local  subdomain,  f22j  while 
the  remaining  unknown  values,  Uj,  are  found  on  the  complement  of  the  subdomain,  Cfi2* 
Applying  a  similar  decomposition  to  the  restricted  relaxation  operator  has  the  form 


)  ’ 


where  is  a  relaxation  scheme  corresponding  to  the  local  problem  on  resulting  in 


changes  to  U2  only.  The  local  multigrid  method  ,  defined  by  can  then 

be  expressed  in  the  two-grid  form  as  follows: 

,f2/.  ^  J2h^fh  _ 

•U2^  -t-  (^2h)-lf2fc 

•  u'*  •<-  u'*  -+- 
•u'*  <-  G^j(u'*,f'*) 

The  objective  of  the  local  multigrid  scheme  MG^^  is  to  eliminate  unnecessary  computations 
on  Gf22-  But  only  the  relaxation  process  has  been  suppressed  there,  not  the  intergrid  trans¬ 
fers  or  the  residual  calculations.  Fortunately,  the  intergrid  transfers  and  residual  calculations 
can  be  suppressed  without  changing  the  final  result.  With  u'*  as  the  initial  approximation 
for  one  cycle,  as  the  correction  computed  in  the  above  two-grid  form,  and  v* 

as  the  change  in  -I-  resulting  from  the  relaxation  on  (again  from  above),  the 

approximation  produced  by  MGq^  is  u^-t-/2/,u^'‘+v'‘.  Since  v/,  is  zero  on  002  (no  relaxation 
occurs  on  CO2  in  the  relaxation  scheme),  is  zero  on  COj.  or  rather  on  the  set  of 

points  that  are  not  in  O2  or  its  interface  (  Oj  0002)-  Using  the  two-grid  form  above  and 
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the  Galerkin  condition  results  in 


^2fc(fh  _  ^  jh^^2h  ^ 


=  /f  (f'*  -  A’'u^)  -  /f  -  ll^A'^v^ 

-  /f  (f'‘  -  A'^u^)  -  -  /f  ylV'* 

_  ^2h  _  £2/i  _ 


-ll^A^ 


v\ 


This  implies  that  the  transferred  residual  in  the  next  cycle  of  MGq^  will  be  zero  at  the 
points  of  Ci72-  In  general,  local  fine  grid  relaxation  has  no  immediate  effect  on  the  coarse 
grid  equations  at  the  points  on  CQ^-  There  is  no  need  to  use  the  fine  grid  in  that  region 
[Ref.  4]. 

2.  Bordered  Multilevel  Method 

To  prevent  unnecessary  work  away  from  the  local  fine  grid,  the  role  of  in  0^2 
must  change  from  representing  error  corrections  to  representing  the  current  approximation. 
A  modification  of  the  local  multigrid  scheme  MGq^  uses  this  to  accumulate  corrections. 
Also,  when  using  this  modification  of  MGq^  ,  the  fine  grid  must  have  enough  data  to  compute 

—  A^u*)  at  the  coarse  grid  interface  points.  This  requires  interpolation  on  fij 
the  two  grid  lines  bordering  Figure  18  shows  the  first  border  and  the  extended  local  fine 
grid,  Q2+>  which  includes  both  borders. 

These  ideas  provide  for  a  new  algorithm,  called  the  bordered  mulilevel  scheme,  that 
avoids  all  fine  grid  work  away  from  the  local  fine  grid  lu  this  new  algorithm,  residuals  are 
computed  on  the  extended  local  fine  grid,  which  includes  the  local  fine  grid,  the  boundary 
points  of  the  local  fine  grid,  aind  the  first  border.  The  residuals  are  transferred  to  the 
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w  first  border  point 

local  fine  grid  point 

Figure  18.  The  local  grid  border  (indicated  by  •),  and  the  extended  local  grid 

^2+- 

coarse  grid  points  lying  under  the  fine  grid  and  the  interface  (residuals  are  also  computed 
on  the  coarse  grid  outside  the  fine  grid).  A  global  coarse-grid  correction  is  computed  and 
interpolated  to  the  extended  local  fine  grid.  Relaxation  is  then  performed  at  fine  grid  interior 
points. 


3.  Multilevel  Composite  Grid  Method 

The  bordered  multilevel  scheme  requires  two  additional  grid  lines  bordering  the  local 
fine  grid  in  order  to  compute  the  residual  correction  at  the  coarse  grid  interface  points.  This 
residual  correction  is  the  composite  grid  residual  at  the  coarse  grid  interface  points.  Let 
Jlh  denote  the  composite  grid,  U  By  generating  a  composite  grid  operator, 

Ah  :  IJh  Uh  where  IJh  is  the  nodal  space  associated  with  the  composite  grid 
residual  correction  can  be  computed  more  directly  without  the  need  for  any  borders.  This 


leads  to  the  multilevel  composite  grid  method ,  represented  by  ■<—  MZA(uh,  and  defined 
as  follows: 

•Set  =  /f  (fh-ylAuh). 

•Compute  u^'*  ^  (yl2/‘)-if2\ 

•Compute  uh  •<—  uh  -I- 
•Set  =  4(fh  -  A^uh). 

•Set  -f-  0  (the  initial  guess  for  the  fine  grid  is  zero)  and  compute  <—  G^{u^,  f^). 
•Compute  u-  •«— 

In  developing  this  multilevel  composite  scheme,  it  has  become  necessary  to  define 
some  new  notation.  First,  the  domain  Q,  needs  to  be  partitioned  as 

=  (IV.15) 

which  breaks  the  domain  into  subdomsdns.  Qc  is  that  part  of  the  domain  that  contains  all 
coarse  grid  points  excluding  the  interface  points,  17/  includes  the  interface  points,  and 
includes  the  fine  grid  points.  On  the  composite  grid  in  Figure  19,  the  grid  17^  is  partitioned 
into  the  following: 

f2^  =  fi^Uf77Uf7|.  (IV.16) 

f7^  consists  of  the  coarse  grid  points  (excluding  the  interface  points) ,  Qj  includes  the  coarse 
grid  interface  points,  and  fTp  is  the  fine  grid  patch. 

A  block  decomposition  of  the  operators  and  grid  equations  results  from  the  partition 
in  (IV.16).  The  block  decomposition  of  the  grid  equations  A-Ur-  =  fh  is  given  by 
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Figure  19.  Coinposite  grid  pavtitiou.  Subgrid  coutuius  only  coarse  grid  points y  subgrid 
nfj  contains  only  interface  points,  and  subgrid  contains  only  fine  grid  points. 


The  notation  signifies  that  the  equation  being  solved  corresponds  to  an  X-grid 
point  and  the  stencil  reaches  to  a  K-grid  point.  If  and  11^  represent  the  grid  transfer 
operators  restricted  to  the  local  fine  grid,  then  certain  observations  result  from  the  way  the 
grid  operators  are  constructed.  For  the  transfer  equation  the  transfer  operator 

is 


I- 


0/0 


(IV.  17) 


VO  *  4  7 

which  may  be  interpreted  as  follows:  outside  the  local  fine  patch,  interpolation  is  simply 
the  identity.  Interpolating  inside  the  local  fine  patch  is  accomplished  using  standard 
prolongation  away  from  the  borders.  The  star,  ★,  represents  actual  storage  of  the  slave 
points  as  values. 
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The  transfer  operator  is  given  by 


j2h  _ 
lu  — 


/  /  0  0  \ 

0  /  *  .  (IV.18) 

VO  0  / 

Residuals  outside  the  local  fine  patch  are  transferred  by  injection  to  correct  the  equations  on 
the  coarse  grid  U  These  residuals  do  not  affect  Qf  equations.  Rather,  residuals  in 
cire  used  to  correct  equations.  Residuals  in  ^  ^ure  used  to  correct  equations. 
The  transfer  operator  is  given  by 


4"  = 


0 

0 

/ 


(IV.19) 


Residuals  are  transferred  from  the  composite  grid  to  the  fine  grid  patch  in  the  trivial  manner. 
The  transfer  operator  is  given  by 


4  =  (»  »  0- 


(IV.20) 


Interpolation  from  the  fine  grid  patch  to  the  composite  grid  is  the  trivial  interpolation,  the 
identity. 

The  composite  grid  operator  is  given,  as  before,  by 


/  aA  aA  aA  \ 
I  ^cc  ^C1  ^CF  ' 

aA  aA  aA 
^IC  ^II  ^IF 

V  ^FC  -^F/  -^FF  ) 


(IV.21) 


Recall  that  the  notation  A^y  means  that  the  equation  being  solved  corresponds  to  an  X-grid 
point  and  the  stencil  reaching  to  a  T-grid  point.  Since  the  coarse  grid  stencils  do  not  reach 
any  fine  grid  points,  Acpc  and  are  zero. 
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The  composite  grid  operator  can  then  be  written  as 


^  -^cc  *  0 

A-  —  •  ★  *  ★ 

\  0  * 


(IV.22) 


The  composite  grid  operator  is  block  tridiagonal  and  agrees  with  the  coarse  grid  operator  on 
Clc  3^d  with  the  fine  grid  operator  on  flp-  The  stars,  ★,  are  non- zero  but  unspecified  entries 
whose  values  are  not  important  at  this  time. 

The  composite  grid  residuals  in  fic  and  Q,f  can  be  computed  using  the  respective 
coarse  and  fine  grid  stencils.  Special  evaluation  is,  however,  required  at  the  coarse  grid 
interface  points.  The  stencils  for  A-  in  the  multilevel  composite  grid  scheme  MZA  must  be 
computed  for  use  on  fij . 

The  multilevel  composite  grid  scheme  MlA  given  above  is  written  in  the  immediate 
correction  form.  This  scheme  is  considered  to  be  “immediate  correction”  because  both 
intermediate  quantities  and  attempt  to  approximate  the  current  algebraic  error  = 
u^'  —  u-  (u-*  is  the  exact  solution)  while  solving  the  residual  equation 


(IV.23) 


This  immediate  correction  form  attempts  to  solve  the  residual  equation  (IV.23)  by  trans¬ 
ferring  the  composite  grid  residual  to  each  level  in  turn,  computing  a  correction  there,  and 
immediately  interpolating  it  back  to  the  composite  grid. 


4.  The  Fast  Adaptive  Composite  Method 

A  final  modification  to  the  multilevel  composite  grid  scheme  MlAis  made  by  replacing 
fine  grid  relaxation  with  a  direct  solver.  The  replacement  of  the  fine  grid  relaxation  with  a 
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direct  solver  is  given  by 


This  modification  gives  the  fast  adaptive  composite  grid  method  (FAC).  This  method  is 
represented  by  u-  ■<—  FAC-{u-,  {-)  and  is  defined  as  follows: 

•Set  =  /f  (fh  -  A^uA). 

•Compute  <— 

•Compute  ^ 

•Set  f'*  =  /^(f^  -  A^u^). 

•Compute  u*  <- 
•Compute  uh  •<—  uh  +  /^u^. 

C.  DISCRETIZATION  ON  THE  COMPOSITE  GRID 

To  apply  FAC,  the  equations  to  be  solved  must  be  discretized  on  a  composite  grid. 
FVE  provides  the  method  for  this  discretization  since  FVE  has  proven  effective  with  PDEs 
on  a  composite  grid. 

1.  Discretization  of  the  Model  Problem  on  a  Composite 
Grid 

As  an  example,  recall  that  the  model  problem  is  the  potential  equation 

V  •  {pVip)  =  T]  in  n  (IV.24) 

Tp  =  xpo  on  (Dirichlet) 

{pVip)-h  =  '01  on  (Neumann). 

The  composite  grid  can  be  taken  as  a  coarse  grid  with  a  patch  of  fine  grid  points,  where  the 
coarse  grid  points  and  the  slave  interface  points  along  the  boundary  of  the  fine  grid  patch 
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appear  in  the  patch  equations  as  if  they  were  Dirichlet  boundary  points  for  the  patch  [Ref. 
4].  The  patch  volumes  are  chosen  first,  with  the  volumes  along  the  interface  being  largely 
dictated  by  the  neighboring  patch  volumes.  Triangulation  follows  with  the  patch  geometry 
being  determined  at  the  interface  by  combining  pairs  of  triangles  of  opposing  orientation 
that  abut  the  slave  points  (Figure  20). 


Figure  20.  Example  of  composite  grid  volumes  and  composite  grid  triangulation. 
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Figure  21.  Interface  point  volume  for  a  composite  grid. 


As  an  example  of  a  stencil  for  a  composite  grid,  consider  the  interface  point  P  in 
Figure  21.  The  control  volume  is  smaller  for  the  interface  point  than  a  control  volume  for 
the  same  point  on  a  uniform  coarse  grid,  and  there  are  now  seven  triangular  regions  over 
which  to  integrate.  In  the  seven  regions  are  ten  portions  of  the  border  of  the  volume.  Using 
the  same  procedures  as  for  the  general  interior  grid  point  of  the  uniform  grid,  integrate  over 
each  of  the  ten  segments  of  the  volume  boundary.  The  resulting  stencil  form  of  the  equation 
for  the  interface  point  P  is 

^*U+2 

PiJ—^  2Pi+\j—\ 

Here,  as  before,  the  notation  E  represents  the  sum  of  the  surrounding  stencil  entries. 

2.  Discretization  of  the  Weldpool  Problem  on  a  Composite 
Grid 

With  an  understanding  of  how  FVE  can  be  used  to  discretize  a  PDE  on  a  composite 
grid,  now  consider  discretization  of  the  weldpool  problem  on  a  composite  grid.  Equations 
II.5  -  II.8  are  considered  now  on  a  composite  grid.  Figure  22  shows  a  two-level  single  patch 
composite  grid,  with  grid  spacing  on  the  fine  patch  taken  as  half  the  grid  spacing  of  the 
coarse  grid. 

The  composite  grid  points  that  will  be  examined  are  those  along  the  coarse-fine  inter¬ 
face.  Six  grid  points  are  labeled  in  Figure  22.  These  grid  points  me  interior  interface  points 
(P  and  Q),  Neumann  boundary  interface  points  (M  and  N)  cind  Dirichlet  boundary  interface 
points  (D  and  F).  For  the  Dirichlet  boundary,  recall  that  since  the  Dirichlet  conditions  are 


ipij  =  ^T]{ih,jh). 
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Figure  22.  A  simple  composite  grid  with  volumes  around  interior  boundary  points  P  and 
Q,  Neumann  boundary  points  M  and  N,  and  Dirichlet  boundary  points  D  and  F  ( extended 
volumes  around  D1  and  FI). 

imposed  directly  on  T,  the  control  volumes  are  extended  so  that  the  control  volumes  nor¬ 
mally  associated  with  the  boundary  points  are  taken  into  account  by  an  extended  volume 
from  the  neighboring  points  (D1  and  FI). 

Beginning  with  the  conduction  equation  for  the  solid, 

^  (IV.25) 

at 

the  discretization  of  the  conduction  problem  resulted  in  the  following  integrals: 

aMo-'  f  '  •  nd5]  If/'  (IV.26) 

J  s 

t  f;  L  +  (1  -  L  •  "<<5]  TTj. 

ij=0 19  J 
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As  before,  the  integrals  over  the  basis  functions  are  calculated  by  representing  the  basis 
functions  as  collections  of  triangular  planes.  The  triaingular  planes  are  assembled  to  form 
the  “hat”  function.  An  interior  grid  point  along  the  interface  has  seven  triangular  planes 
joined  in  this  manner  (Figure  21).  A  Neumann  boundary  point  has  four  such  triangular 
planes,  as  can  be  seen  in  Figure  23.  A  neighboring  Dirichlet  boundary  point,  with  its 
extended  volume,  also  has  seven  triangular  planes.  (Notice  that  the  edge  labeled  “5'£'3” 
does  not  lie  under  the  “hat”  function  of  the  neighboring  Dirichlet  boundary  point  Dl.) 


Figure  23.  Neumann  boundary  point  (M)  volume  and  Dirichlet  boundary  point  (D)  extended 
volume  on  a  composite  grid. 

Following  the  techniques  outlined  above  for  a  uniform  grid,  interface  point  stencils 
have  been  computed  for  the  volume  and  surface  integrals  of  equation  (IV. 26).  Referring  to 
Figure  22,  the  stencils  have  been  computed  for  interior  interface  grid  points  (P  and  Q)  on  the 
top  and  bottom  boundary,  and  also  for  the  Neumann  (M  and  N)  and  neighboring  Dirichlet 
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V 


(D1  and  FI)  grid  interface  points  on  the  top  and  bottom  boundary.  Below  are  the  stencils 


for  the  volume  and  surface  integrals  for  point  P: 


dy  =  —  32  -  lie 

384 


32 -5c 


192  +  5e 


16  4"  5c 


16  +  6c  Ti 


32  -  10c  48  +  4c  16  +  6c 


•  nj  dS 


16  -  9c  -128 

16  -  6e  32  16  +  6c 


16  +  9c  Ti.. 


In  the  above  stencils,  c  =  ^  (r  is  the  radial  coordinate  for  the  central  point  of  the 
volume),  and,  as  before,  the  27r  resulting  from  the  integration  in  the  0  direction  has  been 
factored  out.  The  stencils  for  the  remaining  boundary  points  can  be  found  in  the  Appendix. 

The  next  stencils  to  be  calculated  are  those  belonging  to  the  convection-diffusion 
equation.  Again,  observe  that  the  discretization  of  the  convection-diffusion  equation, 


-I-  V  •  (uT)  =  Ma-^V^T, 


(IV.27) 


has  the  same  form  as  the  conduction  equation  for  the  solid  except  for  the  term 


[  V-{uT)dV. 

J  V 


(IV.28) 


Interface  point  stencils  have  been  computed  for  the  integral  in  equation  (IV.28).  Referring 
to  Figure  22,  the  stencils  have  been  computed  for  interior  interface  grid  points  (P  and  Q) 
on  the  top  and  bottom  boundary,  and  also  for  the  Neumann  (M  and  N)  and  neighboring 
Dirichlet  (D1  and  FI)  grid  interface  points  on  the  top  and  bottom  boundary.  The  stencil 
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J 


for  point  P  has  the  following  form: 

A  B 

C  D  E 
F  G  H 

where  the  entries  of  the  stencil  are  themselves  stencils  applied  to  the  function  'if.  The  entries 
are  given  here: 


The  stencils  that  have  been  developed  are  for  the  conduction  equation  in  the  solid 
(IV.25)  and  the  convection-diffusion  equation  in  the  liquid  (IV.27).  Also,  only  stencils  for 
upper  interior  interface  grid  points  are  displayed  here.  Many  stencils  have  been  developed 
and  are  displayed  in  the  Appendix. 
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3.  Development  of  the  Weldpool  Problem 

The  stencils  developed  in  the  previous  section  are  for  a  model  of  the  weldpool  problem. 
The  model  used  to  develop  the  above  stencils  assumed  a  stationary  boundary  and  considered 
only  a  flat  interface  between  the  two  states  (Figure  24). 


Figure  24.  An  example  of  a  composite  grid  for  the  simplified  weldpool  problem  with  a  sta¬ 
tionary,  flat  boundary.  The  grid  spacing  is  finer  along  the  solid-liquid  interface. 

The  real  problem  is  far  more  complicated  and  requires  some  method  for  tracking  the 
solid-liquid  interface.  The  next  step  in  solving  the  weldpool  problem  is  to  form  the  boundary 
into  its  naturally  occuring  shape  aind  place  a  fine  grid  over  the  solid-liquid  interface,  where 
the  boundary  is  still  considered  to  be  stationary.  An  example  of  this  composite  grid  is  found 
in  Figure  25. 

The  next  step  is  to  consider  the  movement  of  the  interface.  The  movement  results 
from  the  melting  and  solidifying  of  the  solid.  This  movement  is  governed  by  the  local  heat 
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Figure  25.  An  example  of  a  composite  grid  for  the  weldpool  problem  with  a  stationary,  curved 
boundary.  The  grid  spacing  is  finer  along  the  solid-liquid  interface. 

balance  around  each  interface  point.  A  requirement  for  a  control  volume  along  this  interface 
is  that  both  the  current  time  and  the  next  time  step  must  be  contained  in  the  volume.  Not 
only  the  current  position  of  the  interface,  but  an  approximation  of  the  interface  after  one 
time  step  are  required  to  construct  the  current  local  grid  [Ref.  1].  This  will  require  having 
an  adaptive  component  in  the  time  stepping.  That  is,  at  each  step,  the  next  position  of  the 
interface  must  be  computed.  If  the  change  in  position  is  too  Icirge  and  the  new  interface 
point  lies  outside  the  control  volume,  then  the  time  step  must  be  made  smaller.  A  new  local 
grid  must  also  be  computed  at  each  time  step.  In  this  way,  the  control  volumes  are  chosen 
in  a  manner  that  allows  the  interface  to  be  moved. 

These,  and  other  complicating  factors,  indicate  that  a  good  dead  of  work  remains 
before  the  full  weldpool  problem  can  be  considered  “solved” . 
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V.  CONCLUSION 


The  research  for  which  this  work  is  a  part  is  designed  to  develop  efficient 
numerical  methods  for  solving  the  full  weldpool  problem.  This  work  was  undertaken  with 
the  goals  of  developing  the  FVE  method  for  spatial  discretization  and  finite  differences  for 
time  differentiation  to  discretize  the  equations  that  describe  the  weldpool  problem,  and 
develop  stencil  equations  for  use  in  FAC.  This  task  began  with  the  application  of  finite 
differences  to  the  time  derivatives  in  the  conduction  problem,  resulting  in  a  time-stepping 
scheme  that  makes  use  of  a  weighted  average  at  the  current  and  the  next  time  steps.  FVE  is 
then  applied  to  the  potential  flow  problem  to  develop  the  method  for  discretizing  the  spatial 
derivatives.  Having  developed  these  methods  on  model  problems,  the  weldpool  problem  is 
tackled  by  applying  finite  differences  and  FAC  to  discretize  the  governing  equations. 

Multigrid  is  presented,  including  a  development  of  the  coarse  grid,  intergrid  transfers 
and  coarse  grid  correction.  Multigrid  is  introduced  here  because  of  its  effectiveness  as  an 
iterative  method.  Multigrid  proves  sufficient  for  solving  problems  on  a  uniform  grid.  How¬ 
ever,  the  weldpool  problem  requires  an  extremely  fine  level  of  accuracy  along  the  liquid-soUd 
interface.  To  achieve  this  level  of  accuracy  on  a  global  grid  would  be  too  costly.  Thus,  FAC 
is  developed  to  make  use  of  composite  grids.  These  grids  allow  finer  grid  spacing  along  the 
boundary  of  the  molten  pool  and  the  solid,  while  still  providing  for  coarser  grid  spacing 
where  fine  accuracy  is  not  required.  To  arrive  at  FAC,  several  ideas  are  drawn  from  sim¬ 
pler,  less  efficient  methods  related  to  multigrid,  including  local  multigrid  and  the  multilevel 
composite  grid  method.  A  simple  model  of  the  weldpool  problem  is  then  discretized  on  a 
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composite  grid  using  FVE.  Stencil  equations  are  developed  for  the  conduction  equation  in 
the  solid  and  the  convection-diffusion  equation  in  the  liquid. 

Finite  Volume  Element  methods  are  shown  to  be  an  effective  method  for  discretizing 
equations  on  a  composite  grid.  The  stencil  equations  that  result  from  the  discretization  will 
be  coded  into  FAC  at  a  later  date  to  determine  their  effectiveness.  Future  work  will  center  on 
applying  the  FVE  discretization  and  FAC  to  a  moving,  time-dependent  liquid-solid  interface. 
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APPENDIX.  STENCIL  EQUATIONS 


Included  here  are  the  remaining  stencil  equations  that  were  not  included  in  the  main 
body  of  the  text. 

The  first  of  the  stencil  equations  are  for  the  conduction  problem.  The  discretization 
of  the  conduction  problem  resulted  in  the  following  integrals: 


Figure  26.  A  composite  grid  with  volumes  around  interior  boundary  points  P  and  Q,  Neu¬ 
mann  boundary  points  M  and  N,  and  Dirichlet  boundary  points  D  and  F  ( extended  volumes 
around  D1  and  FI). 


Referring  to  Figure  26,  the  stencils  have  been  computed  for  the  conduction  problem 
above  (.1)  at  interior  interface  grid  points  (P  and  Q)  on  the  top  and  bottom  boundary,  and 
also  for  the  Neumann  (M  and  N)  and  neighboring  Dirichlet  (D1  and  FI)  grid  interface  points 
on  the  top  and  bottom  boundary.  Given  here  are  the  stencils  for  the  volume  and  surface 
integrals  for  points  Q,  Dl,  FI,  M,  and  N. 


At  point  Q: 


L  i'£Tkjh.l\dV  = 

\k,l 


h'^r 

384 


16  -  6c  48  -  4e  32  +  10c 
16  -  6c  192  -  5c  32  +  11c 


16 -5c 


32  + 5c 


‘‘v7- 


At  point  Dl: 


At  point  FI: 


dV  = 


h^r 


32 -5c 

32 -11c  224  + 26c 

32  -  10c  48  +  4c  48  +  30c 


Tij. 


h‘^r 


16  -  6c  48  -  4c  48  +  20c 
16  -  6c  208  +  5c 


16-  5e 


48  +  49c 


'■ij- 


At  point  M: 


h^r 

3M 


8  +  c  16  +  5c 

104  +  24e  16  +  6e 

32  +  6c  16  +  6c 


Tij. 


At  point  N: 


r  (  \  h'^r 

/  {'^Tki(l)ki  dV  =  — 


\k,l 


16  +  2c  32  +  lOe 
88  + 19c  32  + 11c 

24  + 6c 


T  ■ 

-‘ij- 


72 


At  point  Q: 


/  (Er«V4,r 

•'^'•3  \k,l 


n]  dS  =  —  16  —  9c 

/  32 


16  -  6c  32  16  +  6c 

-128 


16  + 9c  Tij. 


At  point  D1  (Notice  here  that  Tij  is  known  for  the  furthest  right  entries  of  this  stencil. 


These  three  stencil  elements  can  be  pulled  from  the  stencil.): 


16  +  12c 


=  —  16  -  9c 


-128 -9c  -32 -18c  Tij. 

16  -  6c  32  16  32  +  30c 


At  point  FI  (Notice  here  that  Tij  is  known  for  the  furthest  right  entries  of  this  stencil. 
These  three  stencil  elements  can  be  pulled  from  the  stencil.): 


16  -  6c  32  16  32  +  30c 


16 -9c 


At  point  M: 


-112  +  7c 


16 -6c 


-16  -  10c  Tij. 


16  +  4c 

-64 -21c  16  + 9c  Tij. 

16  +  2c  16  +  6c 


At  point  N: 


•'^'•3  \k,l 


n  1  d5  = 


16  +  2e  16  +  6c 
-64  -  21c 


16  +  4c 


16  +  9c  Ti 


In  the  above  stencils,  c  =  ^  (r  is  the  radial  coordinate  for  the  central  point  of  the 
volume)  and,  as  before,  the  27r  resulting  from  the  integration  in  the  6  direction  has  been 


factored  out. 
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Discretization  of  the  convection-diffusion  equation, 


— -hV-(uT)  =  Ma-iv2T,  (.2) 

has  the  same  form  as  the  conduction  equation  for  the  solid  except  for  the  term 

^V-(uT)dV.  (.3) 

Interface  point  stencils  have  been  computed  for  the  integral  in  equation  (.3).  Again  referring 
to  Figure  26,  the  stencils  have  been  computed  for  interior  interface  grid  points  (P  and  Q) 
on  the  top  and  bottom  boundary,  and  also  for  the  Neumann  (M  and  N)  and  neighboring 
Dirichlet  (D1  and  FI)  grid  interface  points  on  the  top  and  bottom  boundary.  The  stencils 
given  here  ard  for  points  Q,  Dl,  FI,  M,  and  N. 

At  point  Q,  the  stencil  has  the  following  form: 

ABC 
DBF 

G  H 

where  the  entries  of  the  stencil  matrix  are  themselves  stencils,  applied  to  the  function 
The  stencil  entries  are  given  as: 
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At  point  Dl,  the  stencil  has  the  following  form: 

/  A 


|v.(«T)dV  =  1 


B  C 
D  E  F 


J 


where  the  entries  of  the  stencil  matrix  are  themselves  stencils,  applied  to  the  function 
The  stencil  entries  are  given  as: 
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At  point  FI,  the  stencil  has  the  following  form: 
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where  the  entries  of  the  stencil  matrix  are  themselves  stencils,  applied  to  the  function  'I' 
The  stencil  entries  are  given  as: 


At  point  M,  the  stencil  has  the  following  form: 

(A  B\ 

/,V.(uT)dV  =  1  ?  ,  ^ 

h  r 

V  y 

where  the  entries  of  the  stencil  matrix  are  themselves  stencils,  applied  to  the  function 
The  stencil  entries  are  given  as: 
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At  point  N,  the  stencil  has  the  following  form; 


A  B 


y^V-(uT)dV  =  —  CD  Tij, 


where  the  entries  of  the  stencil  matrix  are  themselves  stencils,  applied  to  the  function 
The  stencil  entries  are  given  as: 
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The  stencils  that  have  been  included  here  are  for  the  conduction  equation  in  the  solid 
(IV.25)  and  the  convection-diffusion  equation  in  the  liquid  (IV.27).  The  stencils  here  are  for 
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the  interface,  either  in  the  interior  or  along  the  boundary.  Stencil  equations  for  the  vorticity 
and  stream  functions,  equations  (II. 7)  and  (II.8),  still  remains  for  future  computation. 
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